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Preface 


This  study  was  part  of  a  continuing  effort  to  design  a 
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laser  weapons  using  modem  estimation  techniques.  An  adap¬ 
tive  extended  Kalman  Filter  which  uses  direct  infrared 
sensor  data  was  synthesized  and  tested.  We  hope  that  this 
effort  brings  AFWL  closer  to  an  online  realization  of  such 
a  filter. 
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Abstract 


An  adaptive  Extended  Kalman  Filter  algorithm  is  designed 
to  track  a  distributed  (elliptical)  source  target  in  a  closed 
loop  tracking  problem,  using  outputs  from  a  forward  looking 
infrared  (FLIR)  sensor  as  measurements.  The  filter  adap¬ 
tively  estimates  image  intensity,  target  size  and  shape, 
dynamic  driving  noise,  and  translational  position  changes 
due  to  two  effects:  actual  target  motion,  and  atmospheric 
jitter.  Atmospheric  backgrounds  are  studied  for  the  effect 
of  temporal  and  spatial  correlations  on  filter  performance. 

A  Monte  Carlo  analysis  is  conducted  to  determine  filter 
performance  for  two  target  scenarios:  approximately  straight 
approach  and  cross  range  constant  velocity.  Good  perform¬ 
ance  is  obtained  for  the  first  two  trajectories.  For  the 
second  trajectory,  a  one  sigma  tracking  error  of  .2  pixel 
(4  yrad)  with  a  signal  to  noise  ratio  of  12.5.  The  filter 
adapts  well  to  changes  in  image  intensity,  size,  and  shape. 


AN  ADAPTIVE  DISTRIBUTED-MEASUREMENT 
EXTENDED  KALMAN  FILTER  FOR  USE  IN  A 
SHORT  RANGE  TRACKER 


I.  Introduction 

"...I  believe  that  if  the  difficult  technical . hurdles 
facing  directed  energy  technology  are  surmounted,  the 
application  of  this  technology  to  military  needs  may  revo¬ 
lutionize  both  strategic  and  tactical  warfare."  (Ref  1:48) 

Dr.  Ruth  M.  Davis,  Deputy  Under 
Secretary  of  Defense  for 
Research  and  Advanced  Technology 


Background 

Since  its  conception  and  early  development,  the  laser 
has  been  desired  as  a  weapon  because  of  its  unique  charac¬ 
teristics.  Chief  among  these  is  the  speed  of  light  at 
which  energy  is  transmitted  from  source  to  destination. 
Thus,  a  laser  virtually  eliminates  time  of  flight  consider¬ 
ations;  it  is  no  longer  necessary  to  "lead"  the  target  as 
other  weapon  systems,  such  as  radar  guided  anti-aircraft 
artillery,  must.  With  a  laser,  the  destructive  energy 
reaches  the  target  essentially  instantaneously.  Recent 
technological  advances  with  high  energy  lasers  have  made 
the  deposition  of  large  amounts  of  electromagnetic  energy 
on  a  small  area  in  short  periods  of  time  feasible.  This 
increase  in  laser  technology,  combined  with  speed  of  light 
transmission,  has  attracted  the  attention  of  high  level 
military  and  Department  of  Defense  leaders,  such  as  William 
J.  Perry,  Under  Secretary  of  Defense  for  Research  and  Engi¬ 
neering  (Ref  2:62),  and  General  Alton  D.  Slay,  Commander  of 
the  Air  Force  Systems  Command  (Ref  3:50),  as  well  as  Dr. 
Ruth  B.  Davis,  Deputy  Under  Secretary  of  Defense  for 
Research  and  Advanced  Technology,  whose  quote  begins  this 
paper. 
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Implementation  of  this  desirable  weapon  system  requires 
high  power  for  laser  pumping,  very  accurate  optics  control 
(pointing) ,  and  very  accurate  target  position  estimation. 

The  last  two  requirements  are  dictated  by  the  need  to  main¬ 
tain  the  laser  beam  on  a  specific  part  of  the  target  long 
enough  to  destroy  the  vehicle  or  one  of  its  critical  com¬ 
ponents.  "Painting"  the  entire  target  with  laser  energy  is 
inefficient  and  would  require  extremely  high  amounts  of 
energy  to  achieve  destruction,  and  a  very  accurate  pointing 
and  tracking  system  makes  it  unnecessary.  The  focus  of 
this  research  is  to  investigate  a  precise  target  estimation 
algorithm  for  possible  use  in  laser  pointing  and  tracking. 

An  adaptive  extended  Kalman  filter  as  a  missile  tracker 
has  considerable  potential  to  use  more  information  from  the 
sensor  measurements  than  a  correlation  algorithm  can. 

First,  the  statistical  characteristics  of  the  atmospheric 
disturbance  of  the  infrared  (IR)  wavefront  are  well  known 
and  can  be  modelled  in  the  extended  Kalman  filter.  Thus, 
separation  of  the  atmospheric  jitter  from  true  target  motion 
is  possible.  This  separation  is  important  since  the  high 
energy  laser  wavefront  will  not  undergo  the  same  atmospheric 
distortion  as  the  infrared  wavefronts  (Ref  4:1).  Second, 
the  extended  Kalman  filter  has  the  ability  to  predict  future 
target  motion  through  an  internal  target  dynamics  model.  An 
adaptive  extended  Kalman  filter  can  also  provide  fast  detec¬ 
tion  of  target  maneuvering  and  can  estimate  uncertain  param¬ 
eters  such  as  descriptors  of  target  size  and  shape.  Finally, 
as  needed,  the  filter  can  model  the  temporal  and  spatial 
correlations  of  background  noise,  which  would  enhance  the 
separation  of  atmospheric  jitter  and  target  dynamics. 

Problem 

The  Air  Force  Weapons  Laboratory  (AFWL)  located  at 
Kirtland  Air  Force  Base,  New  Mexico,  is  currently  using  a 
Forward  Looking  Infrared  (FLIR)  Sensor  in  conjunction  with 


2 


/ 


a  "coarse"  correlation  tracking  algorithm  and  a  fine  track 
Far  Field  Intensity  Measurement  (FFIM)  system  to  provide 
accurate  target  position  estimates.  Specifically/  AFWL  is 
interested  in  tracking  air-to-air  missiles  at  close  ranges 
to  accuracies  better  than  one  milliradian  standard  devia¬ 
tion.  As  an  alternative  to  the  correlation  tracker,  an 
adaptive  extended  Kalman  filter  is  developed  in  this 
research.  This  Kalman  filter  does  not  process  correlation 
tracker  outputs  but  uses  FLIR  outputs  directly. 

Research  Outline 

In  1978,  Captain  Daniel  Mercier  completed  an  initial 
thesis  which  demonstrated  the  feasibility  of  using  an 
extended  Kalman  filter  that  uses  direct  FLIR  data  for  target 
position  estimation  (Ref  5) .  The  target  he  studied  was  at 
long  ranges,  and  could  be  modelled  as  a  point  source  of 
infrared  (IR)  radiation.  Further,  the  target  could  be  well 
described  dynamically  by  a  zero-mean  first  order  Gauss- 
Markov  process  for  position  in  each  direction  of  the  two 
planar  axes  in  the  FLIR  field  of  view.  The  target  motion 
model  had  a  small  RMS  (root  mean  square)  value  resulting  in 
small  excursions  for  each  sample  period  of  the  FLIR  sensor. 

Mercier' s  analysis  consisted  of  a  "truth  model"  to  represent 
real  world  target  motion  and  disturbance  environment,  and  a 
filter  (tracker)  to  estimate  the  target  motion  from  meas¬ 
urements  generated  from  the  truth  model  (see  Appendix  A  for 
details) . 

Mercier' s  truth  model  and  filter  were  the  starting 
point  for  the  research  presented  in  this  paper.  Initially, 
the  robustness  of  Mercier 's  filter  was  tested  by  running 
computer  simulations  with  truth  model  and  filter  parameter 
mismatches.  Next,  major  truth  model  changes  were  made  to 
accommodate  the  new  close-range,  air-to-air  missile  target 
scenario,  and  Mercier' s  filter  was  again  tested  against  the 
new  truth  model  in  computer  simulations.  The  results  of 
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the  two  sets  of  simulations  indicated  the  areas  of  appro¬ 
priate  study  for  filter  improvement.  This  list  of  improve¬ 
ment  areas  was  prioritized  according  to  their  effect  on 
filter  performance.  Subsequently,  filter  changes  were  made 
in  a  logical  sequence,  and  computer  simulations  were  run  to 
evaluate  new  performance  capabilities.  The  result  is  a 
filter  which  could  potentially  meet  all  AFWL  requirements. 

Having  described  the  research  goals  in  general  terms, 
certain  assumptions  are  required  to  further  define  and  limit 
the  problem. 


Assumptions 

FLIR.  The  FL1R  outputs  are  instantaneous  samples  of 
an  array  of  infrared  detectors  which  are  mechanically 
scanned  by  a  many  facetted  mirror  through  a  restricted  field 
of  view.  Each  minute  infrared  detector  emits  electrical 
current  proportional  to  the  intensity  of  the  infrared  pho¬ 
tons  entering  the  face  of  the  detector.  A  single  digitized 
output  represents  a  real  time,  electronic  spatial  average 
of  n  horizontally  scanned  detectors.  The  serial  digitized 
data  can  either  be  stored  or  displayed  on  a  cathode  ray 
tube  (CRT) ;  each  picture  element  of  the  CRT  is  called  a 
pixel.  The  horizontal  scanning  of  the  detectors  proceeds 
vertically  through  the  FLIR  field  of  view  resulting  in  an 
array  (frame)  of  pixels  which  is  analogous  in  size  (about 
500  by  400  elements)  and  appearance  to  a  normal  TV  picture. 

A  new  frame  of  pixels  is  generated  every  thirtieth  of  a 
second  (30  Hz  frame  rate) .  For  this  study,  an  8-by-8  array 
of  pixels  out  of  each  large  frame  of  pixels  constitutes  a 
single  measurement  array  (often  called  a  tracking  window) . 
Restricting  the  size  of  the  measurement  array  is  primarily 
dictated  by  computational  and  storage  limitations  and  par¬ 
tially  justified  by  fast  measurement  rates  (to  be  analyzed 
later).  (Ref  5:4-5) 


The  digitized  FLIR  data  has  notable  biases  in  the 
alternating  rows  of  data.  The  elements  of  two  rows  (called 
rows  A  &  B)  of  pixels  are  generated  alternately  during  each 
horizontal  scan.  With  a  clear  field  of  view  (i.e.  FLIR 
pointed  at  blue  sky)  the  average  of  an  A  and  B  row  of  out¬ 
puts  will  differ  by  one  or  two  step  counts.  This  quirk  in 
the  data  can  be  easily  compensated  so  it  is  not  modelled  in 
the  study  simulation  program.  However,  this  data  discrep¬ 
ancy  must  be  accommodated  when  analyzing  real  digital  data 
for  spatial  and  temporal  correlations. 

Target.  The  target  for  current  AFWL  tests  and  this 
research  is  an  air-to-air  missile.  At  close  ranges,  a 
missile  will  appear  to  have  its  own  distinct  geometrical 
shape.  Previously,  Mercier  modelled  the  long  range  target 
as  a  point  source,  yielding  a  bivariate  Gaussian  intensity 
profile  with  circular,  egual  intensity  contours  on  the  FLIR 
image  plane.  For  close  range  missiles,  an  extension  to 
this  model  to  allow  the  bivariate  Gaussian  model  to  have 
elliptical,  equal  intensity  contours  with  specified  angular 
orientation  is  used.  Real  FLIR  data  supports  the  use  of 
this  model  to  approximate  closely  missile  shape  in  the  image 
plane  (FLIR  focal  plane) . 

Air-to-air  missiles  can  exhibit  high  velocities  as 
well  as  high-g  turns  -  a  very  different  target  motion  from 
that  which  Mercier  studied.  This  will  require  substantial 
alterations  to  the  assumed  model  of  target  dynamics. 

To  track  an  air-to-air  missile,  its  image  roust  be  kept 
within  the  8-by-8  field  of  view.  For  a  benign  target  (very 
slow  moving) ,  an  estimate  of  only  position  offset  variables 
would  be  needed  to  provide  to  the  controller  as  the  intens¬ 
ity  pattern  moves  less  them  half  the  field  of  view  in  one 
sample  period.  However,  the  speed  of  a  missile  enables  the 
intensity  pattern  to  leave  the  field  of  view  in  one  sample 
period.  As  a  minimum,  target  velocity  is  needed  in  addition 
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to  current  position  offset  to  predict  subsequent  target 
position. 

Intensity  Pattern.  The  intensity  pattern  of  the  mis¬ 
sile's  projection  on  the  FLIR  focal  plane  is  well  modelled 
by  elliptical,  equal  intensity  contours  (Ref  14:15).  These 
elliptical  contours  can  be  described  mathematically  via  a 
bivariate  Gaussian  intensity  function. 

For  this  model  to  be  useful,  several  parameters  must 
be  known:  the  size  of  the  major  and  minor  axes  and  the 
orientation  of  principal  axes  in  the  image  plane.  Figure  1 
illustrates  an  intensity  function  in  FLIR  coordinates. 


I(x,y)  y 


Figure  1.  Gaussian  Bivariate  Distribution  Intensity 
Model,  I(x,y) 


The  equation  for  this  intensity  function  is 

2 


Il*'y)  '  1max  eItp 


-1 

2<l-r)2  . 


(x-’W)  .  ‘y-S'pe.k1 


(x-itPeak>  ^peak1 

a  a 
x  y 


]| 


(1) 


( 


6 


where  Imax  =  the  maximum  target  intensity 

x,y  =  the  coordinates  of  any  point  in  the  FLIR 
focal  plane 

xpeak'ypeak  =  coorcJinates  locating  peak  intensity, 
*max 

r  =  Image  dispersion  correlation 

axiay  -  Image  dispersions  in  the  x  and  y  direc¬ 
tions  respectively 

The  intensity  function  would  appear  as  elliptical,  equal 
intensity  contours  in  the  FLIR  as  shown  in  Figure  2. 


Figure  2.  Target  Image  in  FLIR  Field  of  View 

The  apparent  target  intensity  function  location  in  the 
FLIR  field  of  view  consists  of  several  components.  Bore- 
sight  error,  FLIR  system  vibrations  and  others  are  assumed 
to  be  negligible  so  that  the  intensity  function  location 
can  be  centered  by: 


where  =  position  offsets  due  to  target  dynamics 

xA'yA  =  Position  offsets  due  to  atmospheric  jitter 

The  separation  of  Xp  from  xA  and  yD  from  yA  is  desirable 
for  the  target  tracker,  a  capability  not  possible  with  cor¬ 
relation  trackers. 

Background  Noise.  The  image  background  noise  which, 
along  with  FLIR.  sensor  noise,  contaminate  the  FLIR  measure¬ 
ments,  is  modelled  as  a  spatially  and  temporally  correlated, 
Gaussian  process.  Various  physical  backgrounds  result  in 
FLIR  images  with  differing  spatial  and  temporal  correlation 
characteristic^.  Real  data  analysis  and  AFWL  experience 
has  been  used  to  determine  appropriate  correlation  coeffi¬ 
cients  for  the  spatial  and  temporal  correlations. 

Closed  Loop-.  This  is  a  closed  loop  tracking  system. 

The  laser  poin^.ng  system  is  assumed  to  be  perfect.  That 
is,  the  system  ca*  .goint  exactly  where  the  tracker  commands 
it  within  each  sample  time.  This  implies  that  settling  time 
of  the  pointing  system  is  less  than  the  data  sample  period, 
the  time  between  discrete  samples  of  the  FLIR  output. 


Previous  Investigation 


The  truth  model  and  filter  derived  by  Mercier  will 
serve  as  a  beginning  for  this  research.  The  truth  model 
target  dynamics  were  modelled  as  independent,  first  order 
Gauss-Markov  processes  in  the  two  FLIR  component  directions, 
written  mathematically  ass 


xD(t)  -  (-  l/xD)xD(t)  +  wx(t)  (3) 
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(4) 


yD(t)  =  (-  l/tD)yD(t)  +  w2(t) 

where  w1  and  w2  are  independent  white  Gaussian  noises  and 
td  is  the  dynamics  correlation  time  (refer  to  Appendix  A 
for  a  more  detailed  description  of  Mercier's  work)  (Ref  5s 
9-10) . 

The  truth  model  atmospheric  disturbances  were  modelled 
as  third  order  Gauss-Markov  processes  based  on  studies  by 
The  Analytical  Sciences  Corporation  (TASC)  (Ref  6,7).  As 
part  of  this  study,  TASC  generated  a  model  for  xA,  the 
atmospheric  jitter,  in  equation  (2)  as  the  output  of  a  sys¬ 
tem  with  a  frequency  domain  transfer  function  of 

*A  =  K(14.14) (659. 5)2  (5) 

W3  (s  +  14.14) (s  +  659. 5)2 

as  driven  by  a  unit  strength  white  Gaussian  noise  w3.  K  in 
equation  (5)  is  adjusted  to  establish  the  current  root  mean 
square  (RMS)  value  of  xA.  yA  is  similarly  modelled.  (See 
Appendix  A  for  conversion  of  equation  (5)  to  Jordan  canoni¬ 
cal  form. ) 

These  models  resulted  in  an  eight  state  truth  model 
with  one  dynamic  state  and  three  atmospheric  states  for 
each  direction.  The  truth  model  state  vector  is  propagated, 
after  conversion  to  equivalent  discrete  time  form,  by  the 
following  equations 

2T(ti+l)  =  £t  W  +  W  (6) 

where  S-p^t^),  ^(ti^  =  state  vector  at  times  ti+1  and  t^ 

respectively 

w^(t^)  *  discrete  noise  term 

$T  =  constant  state  transition  matrix  (Ref  5:15) 

The  truth  model  target  intensity  function  is  converted  into 
measurements  by  the  following  relationship: 


9 


Zkl(ti)  A  ff  Xmax  expj"  2  [(x“xpeak  (tiO 2 
p  i-jth  '  g  L 

pixel 

+  (y-ypeak^iO2])  axdy  +  nkl(ti>  +  £kl(ti> 

(7) 

where  =  output  of  the  i-jth  pixel  at  time 

Ap  =  area  of  a  pixel 

x,y  =  coordinates  of  any  point  in  the  pixel 

Xpeak^i^'  =  coordinates  of  the  center  of  the  apparent 
ypeak(ti)  intensity  distribution  with  respect  to  the 
center  of  the  8-by-8  FLIR  field  of  view 
at  time  t^ 

Xmax  =  maximum  target  intensity 

0^  =  dispersion  of  the  Gaussian  intensity 
function 

nJcl(ti)  =  background  noise  term  for  the  i-jth  pixel 

fkl(t^)  =  FLIR  noise  term  for  the  i-jth  pixel 
(Ref  4:17) 

The  filter  derived  by  Mercier  had  four  states,  one 
state  for  target  dynamics  and  one  state  for  atmospheric 
jitter  in  each  direction.  Both  target  dynamics  and  atmos¬ 
pheric  jitter  were  modelled  as  first  order  Gauss-Markov 
processes  driven  by  white  Gaussian  noise  as  follows 


=  -  ^  *D(t)  +  wx{t) 

(8) 

«»— s 

ft 

77  xA(t)  +  w2(t) 

(9) 
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where  td  =  correlation  time  of  target  dynamics 

w1(t),  w2(t)  =  independent  white  Gaussian  noise 
processes 

ta  =  correlation  time  of  atmospheric  jitter,  1/14.14 

Analogous  equations  were  used  for  yD(t)  and  yA(t)  (Ref  5:19) 
Equation  (9)  is  a  reasonable  approximation  of  equation  (5) 
in  view  of  the  discrepancy  of  the  break  frequencies  and  the 
greater  importance  of  the  lower  frequency  characteristics. 
After  deriving  the  state  transition  matrix,  j>p(t,t^),  and 
Qpd(ti),  the  discretized  noise  strength  matrix,  the  filter 
states  and  covariance  matrix  are  propagated  using  standard 
extended  Kalman  filter  forms  (see  Appendix  A) .  Mercier 
used  the  following  initial  conditions  for  his  filter: 


W '  ' 


<  E(t0! 
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where  P (t  ) 
—  o 


o 
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(10) 


initial  covariance  estimate 
initial  estimate  of  dynamics  variance 

initial  estimate  of  variance  of  atmospheric 
j itter 


The  measurement  update  of  Mercier 's  extended  Kalman 
filter  is  more  unique  and  noteworthy.  First  of  all,  for 
the  design  of  the  filter  the  measurements  (64  pixel  outputs) 
are  assumed  to  be  a  nonlinear  function  of  the  filter  states 
plus  noise  as  follows: 

£(ti)  =  hUpO^),  t±)  +  v(tA)  (11) 

where  h  represents  the  effect  of  the  intensity  spread  func¬ 
tion  given  in  (1) .  This  equation  results  in  the  following 
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measurement  update  equations  in  standard  extended  Kalman 
filter  form: 

-‘V  •  Pit^  )  HT(tiHH(ti)  fltL~)  r'ctj)  +  R<ti) ]-1 

(12) 

s(ti+i  *  ilV1  -hlxltf),  tj)]  (13) 

£(ti+)  *  -  K<ti)  H(ti)  P(ti-)  (14) 


where 


H(tA) 


3h(x,t) 

3x 


x 


x(t.“) 


(15) 


However,  equation  (12)  would  require  the  inversion  of  a  64- 
by-64  matrix.  To  avoid  this,  the  inverse  covariance  form 
of  the  update  cycle  was  used  (Ref  8:238-241)  which  requires 
only  two  4-by-4  matrix  inversions: 

P(ti+)  «  lHT(ti)  R”1^.)  H(tA)  +  P"1^"))"1  (16) 


X(t.+)  =  i(ti")+P(ti+)HT(ti)R-1(ti)[z(t.)-h(£(ti’'),ti)) 


(17) 


Truth  Model  Changes 

Capt.  Mercier's  truth  model  was  modified  in  several 
ways  to  accommodate  this  research.  One  change  was  incor¬ 
porated  to  simulate  an  air-to-air  missile  flight  path. 

This  flight  path  was  programmed  by  adding  a  deterministic 
input  to  the  state  equation. 

x(ti+l>  “  * +  Bu(tA)  (18) 

where  B  is  a  constant  input  matrix  and  u(tA)  is  an  equiva¬ 
lent  discrete  time  input  function  depicting  the  trajectory. 
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To  simulate  the  measurements  taken  from  the  FLIR  image 
plane,  the  truth  model  was  modified  to  use  position  and 
velocity  of  the  missile  to  orient  an  ellipsoid  in  3D  space. 
Assuming  the  missile  velocity  vector  is  colinear  with  the 
major  axis  of  the  ellipsoid,  the  major  axis  of  the  ellipsoid 
is  projected  onto  the  FLIR  image  plane.  This  projection  and 
the  minor  axis  of  the  ellipsoid  are  then  used  as  semi-major 
and  semi-minor  axes  of  the  elliptical  equal  intensity  con¬ 
tours  . 

Mercier  modelled  the  disturbances  in  the  measurements 
due  to  the  background  and  FLIR  (represented  by  n^(t)  and 
(t)  in  equation  (7))  as  spatially  uncorrelated,  white 
Gaussian  noises.  However,  most  FLIR  measurements  exhibit 
both  spatial  and  temporal  relationships  depending  on  the 
background  being  observed.  The  spatial  background  is  easily 
handled  by  generating  spatially  correlated  noise  using  a 
measurement  noise  strength  matrix  with  appropriate  cross 
correlation  terms  (off-diagonal)  filled.  The  modelling  of 
the  temporal  relationship  of  these  noises  is  much  more  com¬ 
plicated.  The  model  chosen  was  to  represent  the  time  his¬ 
tory  of  the  noise  for  each  pixel  as  the  output  of  a  first 
order  lag  system  (yielding  exponentially  time-correlated 
output)  driven  by  spatially  correlated  noise.  Thus,  64  new 
states  were  added  to  the  truth  model  to  accomplish  the  sim¬ 
ulation  of  this  effect. 


II.  Truth  Model  Development 


Introduction 

The  truth  model,  a  mathematical  representation  of  the 
real  world  system  to  be  evaluated,  is  developed  with  those 
characteristics  needed  to  model  the  close  range  air-to-air 
missile  problem  accurately.  Starting  with  the  previous 
truth  model  as  described  in  Chapter  I  and  Appendix  A,  the 
truth  model  simulation  package  is  supplemented  with  several 
large  modelling  changes  needed  to  account  for  the  change  in 
the  target  scenarios.  The  changes  incorporated  are  the 
addition  of  spatial  and  temporal  correlations  to  the  FLIR 
background  noise,  the  portrayal  of  the  target  as  an  intens¬ 
ity  pattern  with  elliptical  constant  intensity  contours, 
the  motion  of  air-to-air  missile  at  close  ranges,  and  the 
effects  of  changing  missile  image  size  and  shape  with  vary¬ 
ing  range  in  a  fixed  field  of  view. 

As  an  alternative  to  the  development  of  a  truth  model, 
the  use  of  an  AFWL-supplied  routine  called  REFERENCING  is 
discussed.  This  routine  is  capable  of  reading  actual  site 
test  data  from  tape.  This  software  package  has  several 
capabilities  and  has  considerable  potential  as  a  future 
analysis  tool. 

Spatially  Correlated  Background  Noise 

The  modelling  of  spatial  correlations  in  the  background 
noise  is  conceptually  easy  but  laborious  in  implementation. 
Before  the  mechanics  for  creating  spatially  correlated  back¬ 
ground  noise  are  presented,  the  need  for  this  modelling  is 
discussed. 

The  existence  of  spatial  correlations  in  the  background 
of  the  FLIR  data  has  been  observed  by  AFWL  personnel  and  can 
be  seen  in  real  data  analysis.  Captain  Lewantowicz  of  the 
AFWL  has  recommended  spatial  correlations  for  distances  up 
to  40  microradians  (yrad)  based  on  his  experience  (Ref  9) . 
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Forty  microradians  is  approximately  equivalent  to  two  pixels 
in  the  scale  used  in  this  research. 

Evaluation  of  real  data,  which  is  explained  in  detail 
in  Appendix  B  along  with  tables  and  results,  further  demon¬ 
strates  that  spatial  correlations  exist  in  the  background 
of  the  FLIR  field  of  view.  The  data  analyzed  was  actual 
digital  FLIR  values  taken  while  tracking  an  approaching 
air-to-air  missile.  The  background  for  the  test  started  as 
clear  sky  and  changed  to  ground  clutter  when  the  missile's 
motion  carried  it  below  the  horizon.  The  data  values  ana¬ 
lyzed  were  the  sum  of  one,  four,  or  sixteen  pixels  (i.e. 
from  individual  pixels,  2-by-2  arrays,  or  4-by-4  arrays)  as 
needed  to  maintain  the  target  image  in  the  limited  field  of 
the  data  matrix.  The  results  of  the  analysis,  which  are 
found  in  Tables  B-I  through  B-IV  of  Appendix  B,  are  sum¬ 
marized  here.  The  correlation  coefficients  between  adjacent 
pixels  (first  neighbor)  in  the  horizontal  direction  aver¬ 
aged  0.325  (the  first  eight  entries  of  Table  B-I)  for  values 
which  represented  a  single  pixel's  output.  The  horizontal 
coefficients  of  values  representing  the  sum  of  four  or 
sixteen  pixels  (the  last  six  entries  of  Table  B-I)  averaged 
0.25  -  slightly  higher  than  expected  because  of  the  large 
distance  which  this  spatial  correlation  represents.  Such 
pixel  summations  were  done  near  the  end  of  the  missile's 
flight  when  the  missile  had  gone  below  the  local  horizon. 
Thus,  the  higher  than  expected  correlations  of  the  multiple 
pixel  values  may  have  been  due  to  a  change  of  physical  back¬ 
ground,  a  background  with  more  structure. 

The  computation  of  vertical,  first  neighbor  coeffi¬ 
cients  was  complicated  by  the  bias  in  the  average  of  alter¬ 
nating  rows  (called  A  and  B)  of  data  as  originally  discussed 
in  Chapter  I.  In  fact,  as  can  be  seen  in  the  first  six 
entries  of  Table  B-III,  the  correlation  coefficients  com¬ 
puted  when  comparing  A  and  B  values  are  primarily  negative. 
When  A  and  B  values  are  added  together  and  then  correlated, 


15 


as  occurs  in  the  last  five  entries  of  Table  B-III,  the 
coefficients  were  positive  and  average  0.355. 

Similarly,  correlation  coefficients  between  a  pixel  and 
the  pixel  two  away  (considering  any  three  pixels  in  a  row, 
the  first  and  third  are  compared)  in  the  vertical  and  hori¬ 
zontal  directions  are  computed  from  real  data  (Appendix  B, 
Tables  B-Il  and  B-IV) .  Horizontally,  the  data  was  unexpec¬ 
tedly  uncorrelated  except  when  four  or  more  pixels  were 
added.  This  may  be  due  in  part  to  the  character  of  the 
background  which  was  clear  sky  when  single  pixel  values 
were  used,  as  in  the  first  eight  entries  of  Table  B-II.  In 
the  vertical  direction  when  four  and  sixteen  pixel  sums, 
which  adds  together  A  and  B  values,  were  correlated,  the 
coefficients  (the  last  five  entries  of  Table  B-IV)  averaged 
0.334. 

Thus,  the  real  data  analysis  of  single  pixels  substan¬ 
tiated  the  existence  of  the  spatial  nature  of  the  background 
of  FLIR  data,  but  the  correlations  were  not  found  to  be  as 
extensive  as  expected  due  probably  to  the  sky  background. 
Other  than  the  last  segment  of  this  test  run,  all  available 
data  sources  had  clear  sky  backgrounds.  Consequently,  the 
spatial  contents  of  other  backgrounds  were  not  studied. 

Another  effect  possibly  contributing  to  the  horizontal 
correlations  was  the  horizontal  sweeping  of  the  detectors. 
This  effect  was  not  noticeable  in  the  data  analysis  because 
the  A  and  B  row  effects  did  not  allow  a  good  comparison  of 
vertical  and  horizontal  coefficients.  More  methodical  data 
analysis  of  known  backgrounds  is  needed  to  establish  the 
horizontal  contribution  of  this  effect  and  to  establish  the 
range  of  reasonable  correlations. 

The  generation  of  spatially  correlated  white  Gaussian 
noises  is  accomplished  by  allowing  nonzero  cross  correla¬ 
tions  between  the  measurement  noises  v^  (t^)  in  the  following 
equation . 
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(19) 


5.(ti>  "  +  v(t±) 

The  components  of  v  (t^) ,  which  is  of  dimension  64,  are  no 
longer  independent  as  in  the  original  truth  model  but  are 
spatially  correlated.  The  noise  covariance  matrix  is  then 


longer  diagonal 

but  has  the  general  form 
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where  R  =  64-by-64  measurement  noise  covariance  matrix 

rkl  =  correlation  coefficient  relating  noise  com¬ 
ponents  v^  and  v1  at  a  given  time  t^ 


v  is  assumed  stationary  so  R  is  constant  for  a  given  back¬ 
ground;  it  would  be  changed  to  reflect  a  significant  change 
in  the  physical  background.  The  correlation  coefficient  is 
defined  statistically  for  zero  mean  random  variables  as: 


k,l 


E^vk^ti^vl^ti) * 

VE{vk(ti)}E{vl(ti)} 


E*vk(ti)vl(ti)  * 

Vi 


(21) 


Prom  the  above  equation,  it  can  be  seen  that  r^ 


°k#  °1» 


V 


and  v^  are  all  scalars  and  that  r^ 


*  rljc  since 
1 .  These 


facts  and  the  fact  that  the  variances  of  the  background 

2  2 

noises  are  assumed  to  have  the  same  strength  (o^  *  c2  “ 
«  a64  “  °l)  yields  a  simplified  noise  covariance  matrix: 
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6 


rl,2 

1 

r2 , 3 

•**  r2,64 

rl,  3 

r2, 3 

1 

r3,64 

rl,64  r2 ,64 


3,64 


(22) 


At  this  point  the  relationship  between  models  for  corre¬ 
lated  and  uncorrelated  spatial  backgrounds  is  most  readily 
understood.  For  uncorrelated  Gaussian  noises,  r^  ^  in  (22) 
equal  zero.  However,  if  the  noises  are  spatially  correlated 
the  noise  in  pixel  k  is  correlated  due  to  its  physical 


proximity  to  pixel  1,  then  some  or  all  of  r^  in  (22)  are 


nonzero. 


Having  explained  what  is  meant  by  spatially  correlated 
noise,  the  next  step  is  to  determine  the  amount  of  spatial 
correlation  that  exists  in  actual  data  and  then  to  reflect 
this  in  R.  AFWL  inputs  and  real  data  analysis  have  shown 
that  the  physical  limit  of  such  correlations  are  one  or  two 
pixels  distance  (first  and  second  neighbors) .  If  the  pixels 
in  the  8-by-8  measurement  array  are  numbered  1  through  64 
starting  in  the  upper  left  and  progressing  left  to  right, 
row  by  row  (see  Figure  3) ,  and  if  the  pixel  noises  are 
similarly  numbered,  then  the  correlation  coefficients  can 
be  assigned  appropriate  values.  Shown  in  Figure  4  are  the 
nonzero  correlation  coefficients  for  pixel  k  for  correla¬ 
tions  limited  to  its  first  neighbor.  The  number  in  the 
upper  right  hand  corner  is  the  pixel  number  with  respect  to 
pixel  k.  The  coefficient  in  each  box  is  that  which  relates 
the  noise  in  that  boxx with  the  noise  in  pixel  k.  For  first  . 
neighbor  correlations  all  coefficients  are  zero  except 
those  relating  pixel  k  to  its  eight  immediate  neighbors. 

The  real  data  analysis  indicated  that  the  coefficients  were 
symmetric  vertically  and  horizontally;  the  four  diagonal 
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Figure  3.  Pixel  Numbering  Scheme 
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Figure  4 .  First  Neighbor  Correlation  Coefficients 


coefficients  were  set  equal  to  each  other,  as  were  the  four 
neighbors  above,  below,  and  to  the  sides.  Thus,  r^ 

rk  ]c-7 '  rk  k+7  ’  and  rk  k+9'  shown  in  Figure  4,  were  set 
equal  to  .2545  for  the  first  neighbor  exponential  curve  fit 

to  be  discussed  subsequently,  and  r^  k_8,  r^  ,  r^  k+1, 
and  rk  }c+g  were  set  equal  to  .3800  for  the  first  neighbor 
exponential  curve  fit.  Making  similar  observations  of  the 
correlation  coefficients  for  all  the  pixels,  the  R  matrix 
can  be  computed. 
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Similarly,  if  the  spatial  correlations  are  to  be 
modelled  as  extending  to  first  and  second  neighbors,  then 
the  coefficients  relating  a  pixel's  noise  to  that  of  its 
24  closest  neighbors  are  nonzero  as  shown  in  Figure  5,  and 
all  others  are  zero.  Making  appropriate  selections  for  all 
pixels  allows  R  to  be  computed,  but  this  time  it  has  con¬ 
siderably  more  nonzero  entries.  Actual  implementation  in 
the  computer  simulation  software  (see  Appendix  G  for  com¬ 
puter  code  listings)  of  the  R  formulation  is  simplified  by 
the  symmetry  of  the  coefficients  about  the  main  diagonal 
and  by  the  use  of  recursive  loops. 
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Figure  5.  First  and  Second  Neighbor  Correlations 


Once  R  is  known,  specific  spatially  correlated  reali 
zations  of  v  can  be  generated  by  using  a  Cholesky  square 


root  decomposition  of  R  and  a  Gaussian  noise  generator  with 
unit  variance.  Specifically,  v(t^)  can  be  generated  as 

vlt^  -  Cy/R  w(t±)  (23) 

where  w(t^)  =  64  dimension  vector  of  independent  white 

noises  of  unit  variance  generated  by  inde¬ 
pendent  calls  to  a  Gaussian  noise  routine, 
i.e.,  E{w(ti)w (t j ) }  =  I  6i j . 

and  *yl[  is  the  Cholesky  square  root  of  R  as  discussed  in 
Chapter  7  of  Reference  10.  The  Cholesky  decomposition  pro¬ 
duces  a  lower  triangular  matrix  such  that 

Vs  Vi?  ■  r  (24) 

The  covariance  of  v  as  given  in  (20)  is  preserved  with  the 
use  of  this  decomposition  as  follows: 

E{v(ti)vT(tj)}  =  E{cV^w(ti)wT(tj)  VlT> 

*  VI  E{Si(ti)wTtt;.)  }  Vl?  =  Vi  I  V^  6ij 

=  R  6ij  (25) 

The  final  step  in  the  spatial  correlation  model  is 
establishing  the  appropriate  values  for  the  correlations. 
The  early  trial  runs  and  real  data  analysis  provided  con¬ 
siderable  insight  into  the  formulation  of  a  good  model  of 
the  correlation  curve.  Initial  runs  were  attempted  with 
high  (r  values  ranging  from  .95  to  .8)  first  and  second 
neighbor  correlations  in  R  and  zero  correlations  elsewhere, 
which  is  physically  an  unrealistic  trend  for  the  spatial 
correlations  to  follow.  This  resulted  in  R  being  numeric¬ 
ally  ill-conditioned. 
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The  Cholesky  decomposition  of  R  could  not  be  computed 
as  R  is  required  to  have  non-negative  eigenvalues  (10:7-5) 
as  computed  on  a  fixed  word  length  computer.  The  Cholesky 
decomposition  will  theoretically  work  for  any  symmetric, 
positive  semidefinite  matrix  given  an  infinite  word  length 
computer.  Because  all  sensors  have  some  noise,  R  should  be 
positive  definite.  Thus,  special  care  must  be  taken  in 
choosing  correlation  coefficients  that  make  sense  physically 
so  that  R  has  numerically  computable  positive  eigenvalues. 
The  real  data  and  the  practical  experience  of  Dr.  P.  S. 
Maybeck  suggested  that  the  correlation  could  be  modelled 
adequately  as  a  decaying  exponential.  The  real  data  analy¬ 
sis  indicated  first  neighbor  correlations  of  about  .3  to  .4 
and  second  neighbor  correlations  of  essentially  zero. 
Remembering  that  autocorrelations  yield  unity  correlation 
coefficients,  these  three  points  (1,0;  .35,1;  .0,2)  suggest 
an  exponential  decay  or  other  nonlinear  curve  fit.  The 
exponential  decay  model  with  correlation  distance  of  one 
pixel  was  chosen  because  it  fit  the  data  for  the  first 
neighbor  and  was  essentially  zero  for  the  second  neighbor 
and  beyond.  The  exponential  model  also  prevented  R  from 

being  ill-conditioned.  There  was  an  upper  limit  found  for 

% 

the  first  neighbor  correlations  and  for  the  first  and  second 
neighbor  correlations,  above  which  R  became  ill-conditioned. 
For  first  neighbor  correlations  this  value  was  about  .4,  and 
for  first  and  second  neighbor  correlations  it  was  about  .5 
for  the  closest  neighbor  and  .3  for  the  second  neighbor. 

In  summary,  the  following  algorithm  can  be  used  to 
generate  a  vector,  v(t^) ,  of  spatially  correlated  noises: 

(1)  Decide  on  the  appropriate  curve  fit  for  the  spatial 
correlations  and  compute  R. 

(2)  Compute  <y/R. 

C3)  Generate  w(t^) ,  the  vector  of  independent  white 
noises  of  unit  variance. 

C4)  Compute  v(t^)  by  premultiplying  w(t^)  by  <ylR. 
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Elliptical  Intensity  Contour  Target 

The  introduction  outlined  in  a  general  fashion  the 
desirability  of  using  elliptical  constant  intensity  contours 
to  model  an  air-to-air  missile  image.  There  are  several 
good  reasons  for  doing  so.  Studying  frames  of  data  from 
AFWL  indicated  that  for  any  particular  intensity  level  the 
pattern  was  approximately  elliptical  and  in  some  cases  cir¬ 
cular.  This  seems  to  bear  out  that  the  atmosphere  and 
optics  tend  to  smear  a  well  defined  image  and  attenuate  the 
edges  gradually.  The  diffraction  pattern  produced  by  the 
image  closely  approaches  a  Gaussian  intensity  function 
(Ref  14:15).  The  use  of  a  bivariate  Gaussian  intensity  func¬ 
tion  to  simulate  this  effect  maintains  the  essence  of  the 
observed  image  properties  and  mathematical  ease  of  handling. 

Capt.  Mercier  assumed  the  object  to  be  far  enough  away 
so  as  to  appear  as  a  point  source  of  radiation.  Due  to  the 
physics  of  wave  propagation  and  optics,  this  yields  an 
intensity  pattern  on  the  FLIR  image  plane  well  approximated 
by  a  bivariate  Gaussian  function  with  circular  equal  intens¬ 
ity  contours.  To  simulate  closer  range  targets  the  point 
source  assumption  was  clearly  inadequate.  This,  along  with 
observed  data,  motivated  the  elliptical  shape  with  speci¬ 
fied  orientation  in  the  image  plane  used  in  this  research. 

The  intensity  function  used  is 


I (x,y) 

« 

-  Vx^P'-T1  <x-xpeak^  ‘^peak’ J  IPI'1  [  (x-xpeak>  ^peak* )T) 

(26) 

where  the  variables  are  shown  in  Figure  6;  xpeak  and  ypeak 
locate  the  peak  of  the  intensity  function,  and  oy  and  0pV 
are  the  eigenvalues  of  the  P  matrix  in  the  equation  above. 
To  quantify  the  image,  a  pixel  dimension  of  20  yrads  is 
used.  Therefore,  an  8  x  8  FOV  will  be  160  yrads  wide  in 
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Figure  6 .  Image  Intensity  Characteristics 

azimuth  and  160  yrads  wide  in  elevation.  Because  of  the 
small  FOV,  angular  displacement  of  the  target  from  the  cen¬ 
ter  of  the  FOV  is  closely  approximated  by  linear  displace¬ 
ment  on  the  FLIR  image  plane.  Also,  angular  velocity 
closely  approximates  image  plane  linear  velocity.  Since 
the  distinction  between  angular  displacement  and  linear 
displacement  is  negligible,  angular  measurements  were  used 
in  this  research. 

The  missile  image  was  produced  in  several  steps  by  the 
truth  model.  From  the  simulated  three-dimensional  inertial 
position  and  velocity  of  the  missile,  the  azimuth  velocity 
(x  velocity  in  Figure  6) ,  elevation  velocity  (y  velocity  in 
Figure  6)  and  speed  (magnitude  of  velocity  vector  in  3-D 
space)  were  computed  in  pixels/sec.  The  azimuth  velocity 
and  elevation  velocity  then  define  the  missile  velocity 
component  perpendicular  to  the  line  of  sight  (LOS) .  The 
ratio  of  the  velocity  perpendicular  to  the  LOS  and  speed  is 
the  cosine  of  the  out  of  plane  angle  y.  Figure  7  shows  the 
geometry  involved. 


With  cosy  determined,  the  semimajor  axis  of  the  three 
dimensional  missile  is  projected  onto  the  image  plane.  This 
projected  length  becomes  ov  and  the  radius  of  the  circular 
missile  IR  cross-section  is  retained  as  the  semiminor  axis 
magnitude,  <j  ,  in  the  image  plane.  With  these  parameters 
established,  the  intensity  at  any  point  in  the  image  plane 
can  be  computed.  This  computation  is  performed  in  the 
principal  axis  coordinate  system  of  the  image  ellipse.  The 
intensity  function  then  becomes 


Ilx'i’1  *  Iroax  exP 


{-  °2  ]  t4xv4yv]T|  (29) 


where  Axy  =  (x-Xpea}?)  cos0  + 

(y-ypeak)sin0 

(30) 

4yv  -  <y-ypeak,cos8  - 

(x-xpeak)sine 

(31) 

where  the  rotation  angle  0  between  (x,y)  and  principal  axis 
coordinates  is  given  in  Figures  6  and  7. 


The  average  intensity  for  any  pixel,  as  measured  by  the 
FLIR,  is  the  integral  of  the  intensity  function  over  the 
pixel  area  divided  by  the  pixel  area.  To  approximate  this 
integral,  the  intensity  was  averaged  over  sixteen  equally 
spaced  points  inside  the  pixel  (i.e.  at  the  centers  of  a 
4-by-4  array  of  squares  in  each  pixel) . 


To  complete  the  average  intensity  measurement  simula¬ 
tion,  noise  due  to  various  FLIR  and  background  sources  was 
added  to  each  pixel.  Finally,  the  intensity  function, 
including  noise  terms,  for  one  of  the  64  pixels  is 


2lm(ti} 


4 

j^max 


4tAxvkdyVj3 


[ixvk*yVj 


"lro^i* 
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(32) 


where  the  noise  term  n^m(t^)  can  k©  modelled  several  differ¬ 
ent  ways.  Three  models  were  used  in  this  research. 

1)  Spatially  and  temporally  white  Gaussian  discrete-time 
noise  where 

E Cn  (t^ )  ]  =  0  (33) 

E[n (t^)nT (t^) ]  =  P  (P  a  diagonal  matrix) (34) 

E(n(ti)nT(tj)]  =  0  i  *  j  (35) 

2)  Spatially  correlated,  temporally  white  Gaussian  noise 
where 

EIMt^]  =  0  (36) 

E[n(ti)nT(ti) J  =  P  (P  not  a  diagonal  matrix) (37) 

EIntti)nT(tj)J  -  0  tL  ft  tj  (38) 

3)  Spatially  and  temporally  correlated  Gaussian  noise  where 

EtnOt^J  =  0  (39) 

Eln  (t^)nT(t^)  J  *»  P  (P  not  a  diagonal  matrix)  (40) 

EI*lCti)nT(tj)J  =  P  e’^i^jl^A  C41) 

ta  =*  correlation  time  of  noise 

With  zero  mean  noise  some  pixels  far  from  the  target  could 
have  negative  intensity  values.  This  was  previously  pre¬ 
vented  by  searching  out  the  minimum  value  of  the  intensity 
pattern  and  setting  this  pixel  value  arbitrarily  to  .1.  The 
other  pixels  were  adjusted  by  the  same  amount  to  maintain 
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the  correct  image  (Ref  5:90).  This  scaling  was  retained 
for  this  research.  The  image  produced  by  the  collection  of 
64  z^m(t^)  is  presented  to  the  tracker  as  the  measured  data. 


Target  Dynamics  Model 

The  continuous  time  dynamic  model  for  the  air-to-air 
missile  used  in  this  truth  model  simulation  was 


with 


x^  *  a  (t) 

x2  **  6  (t) 

u(t)  =  [a(t)  6(t))T 


(42) 

(43) 


In  vector  form,  this  equation  becomes 

>^(t)  «=  u(t) 


(44) 


where  u(t)  is  a  specified  azimuth  and  elevation  velocity 
time  function  used  as  the  input  to  program  a  particular 
trajectory.  Equation  (44)  has  the  general  solution 


*DCt) 


$D(t,to)x(to) 


*/ 


£0(t,T)BD(T)u(T)dT 


(45) 


where  t  =  initial  time 
o 

tjjtt,!)  =  state  transition  matrix  for  missile 
dynamics 

BgCf)  *  input  matrix  for  dynamics 
Ujj(t)  =  input  function  for  dynamics 


The  atmospheric  disturbance  model  as  described  in  Ref  5 
is  a  stochastic  equation  of  the  form 

x^lt)  =  Ct)  (t)  +  (^(tjw^t)  (46) 
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where  3^  (t) 

*A<t) 

^(t) 

Vtl 


the  six  atmospheric  noise  states 

atmospheric  plant  matrix 

atmospheric  noise  input  matrix 

vector  of  white  Gaussian  noise  inputs  with 
zero  mean 


and 


Eiw^tjwj^t+T)  } 


^(tm-r) 


Equation  (46)  has  the  solution 

at  t. 
x 


over  a  sample  period  starting 


A- 


T ) G^ ( f ) w^ ( T ) dx  (47) 


t. 

x 


The  integral  term  in  equation  (47)  requires  a  statistical 
characterization.  Let 


j*  i+l 


*A(ti+l,T)GA(T)wA(T)dT  (48) 


and 


El^Ad(ti)]  =  ° 


(49) 


El^d(ti)-Ad(ti)^ 


*A(ti+i,T)GA(T)gA(T)^(T)^(ti+i 


2ftd 


t  )  dx 


(50) 


(51) 
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Note  that  g Ad  is  not  a  function  of  time:  the  atmospheric 
disturbances  are  modelled  as  stationary  processes. 

For  the  truth  model  to  provide  x^(t.),  realizations  of 

A  X 

wAd are  required.  This  can  be  accomplished  by  taking 
the  Cholesky  square  root  of  and  postmultiplying  it  by 
a  vector  of  unit  strength  white  Gaussian  noise.  Then  equa¬ 
tion  (47)  can  be  written  equivalently  as 

VW  -  +  VSZ  W  (=2) 


where 


E[^ACti)]  =  ° 


(53) 

E^(ti>^(tj)3  "  I«ij 

A  more  detailed  development  of  can  be  found  in  Ref  10. 

Augmenting  equations  (45)  and  (52)  yields  the  complete 
truth  model  time  propagation  model  as 


S^i+l* 


*(ti+l 


■B,(t.n  r  o  -i 

L  » .  J  LVSIJ  1 


(54) 


where  x(t^)  =  state  vector  of  dynamic  and  atmospheric 
states 


/fci+l 

* (ti+1,T)B(T)dT 


(55) 


and  u^ft^)  is  a  piecewise  constant  function  (constant 
between  sample  times) 


Haiti.) 


a  (t^  +  Ap) 

6(ti  +  ^-) 


(56) 


30 


where 


At 


ti+l 


”  ti 


(57) 


Equation  (56)  is  evaluated  at  the  interval  midpoint  as  an 
approximation  to 


/bi+l 

4 


(T)dT  s  o(t±  +  ^)At 


*i 


/*  i+1  . 

J  6 


6(T)dT  3?  $ttA  +  ^)At 


t. 

1 


The  matrices  of  equation  (54)  are 


(58) 


(59) 
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(60) 


The  3x3  submatrices  in  are  identical  and  propagate  the 
atmospheric  noise  states  for  both  azimuth  and  elevation. 
See  Appendix  A  and  Ref  5  for  further  details  of  these  sub¬ 
matrices.  of  (54)  is  specified  as 
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See  Appendix  D  for  details  of  the  matrix. 

Mow,  all  that  remains  is  to  define  o(t)  and  B(t).  Mis¬ 
sile  position  and  velocity  equations  were  programmed  in  an 
inertial  xyz  coordinate  system.  This  was  done  to  ease  the 
transition  from  a  desired  trajectory  to  simulation  equations 
representing  that  trajectory.  The  truth  model  then  takes 
these  inertial  quantities  and  computes  azimuth  velocity  and 
elevation  velocity. 

Figure  8  shows  the  geometry  involved  in  the  azimuth 
velocity  equation.  This  view  is  looking  down  on  the  earth 
surface . 


z 


Figure  8.  Azimuth  Geometry 


.  z_(t)xT{t)  -x_(t)z_(t) 

°  (t)  *  - i-2 - i - ^ -  (63) 

Z][(t)2  +  Xx(t)2 

yielding  an  a(t)  in  rads/sec.  The  truth  model  uses  a  pixel 
as  the  unit  of  angular  displacement  rather  than  radians 
(1  pixel  =  20  yrads) .  Therefore,  a(t)  is  converted  to 
pixels/sec  by  dividing  by  20xl0"”6  rads/pixel. 

Elevation  velocity  is  computed  similarly.  Figure  9 
shows  the  geometry  involved. 


y(t) 


Xj.z-j-  plane 


where  p  =  range  =  [XjCt)2  +  yx(t)2  +  z-^t)2]1^2 


rh  **  horizontal  range  in  xx,zx  P^ane 
«=  IXjtt)2  +  zx(t)2J1/2 


Elevation  can  be  computed  using  either  a  sine  or  tangent 
function.  The  tangent  function  was  chosen  as  it  was  a  sim¬ 
pler  equation  to  implement  using  quantities  already  evalu- 


*h(t)y1(t)  -  Yzlt) 


jXjftJXjCt)  +  ZjftJZjft) 


8  (t) 


rh(t> 


P(t) 


(65) 


Again  this  is  rads/sec  and  is  converted  to  pixels/sec  as 
before  in  the  azimuth  velocity  equation. 

With  a(t^  +  4p)  and  $(t^  +  inserted  into  equation 
(56)  the  truth  model  propagation  equation  is  complete. 

The  trajectories  simulated  by  the  above  equations  were 
chosen  to  exercise  and  verify  several  aspects  of  filter  per¬ 
formance.  The  profile  used  to  evaluate  the  filter  of  Ref  5 
is  shown  in  Figure  10.  The  view  is  looking  down  toward  the 
earth  surface. 


Figure  10.  Trajectory  1:  Approaching  Target 


The  velocities  were  kept  small  so  as  not  to  leave  the  FOV 
in  one  sample  period.  The  equations  are 


i(t.)  -  i0e-a(ti> 

(66) 

6  -  0 

(67) 

where  a  is  a  suitable  time  constant.  This  was  a  simple  tra¬ 
jectory  that  the  filter  of  Ref  5f  with  no  prediction  capa¬ 
bility,  could  track. 
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The  first  trajectory  used  to  test  the  filter  developed 
in  this  research  was  a  constant  velocity  cross  range  path. 
Figure  11  shows  this  path. 


Figure  11.  Trajectory  2:  Cross  Range  Target 
The  equations  are 
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This  path  is  relatively  benign  and  was  designed  to  check 
basic  filter  performance. 

The  next  trajectory  was  a  constant  velocity  cross  range 
path  with  a  20g  pull  up  at  closest  approach  to  tracker  posi¬ 
tion.  Figure  12  shows  this  path. 


Figure 

12. 

Trajectory  3:  20g  Turn 

Before  the 

pull  up 

the 

equations  are 

*  -1000  m/sec 

(74) 

• 

-  0 

(75) 

• 

zl 

=  0 

(76) 

xo 

=  4000  m 

(77) 

*o 

=  0 

(78) 

zo 

=  40,000  m 

(79) 

After  pull 

up  initiation  the  velocity  equations  are 

*1 

8 

-1000  cosll.96 ( t— 4 ) ] 

(80) 

• 

- 

1000  sin[1.96  (t-4) ] 

(81) 

• 

ZI 

8 

0  . 

(82) 

This  path  was  intended  to  show  the  ability  to  track  a  highly 
maneuvering  missile. 


Temporally  Correlated  Background  Noise 

The  analysis  of  real  data  for  temporal  noise  correla¬ 
tions  ,  which  is  described  in  detail  in  Appendix  C ,  was 
inconclusive  except  for  establishing  the  existence  of  such 
correlations.  The  primary  problem  in  arriving  at  conclusive 
results  was  the  lack  of  real  data  to  analyze  with  backgrounds 
other  than  clear  sky.  The  results  of  the  work  accomplished 
(Table  C-I)  demonstrate  a  consistent  temporal  correlation 
coefficient  from  one  sample  time  to  the  next  of  about  .364 
for  the  clear  sky.  This  value  for  a  time  difference  of  one 
thirtieth  of  a  second,  while  disappointingly  low,  was  not 
unexpected  considering  the  back-round,  which  had  essentially 
no  infrared  structure  to  be  correlated  in  time. 

In  the  absence  of  guidance  from  the  data  analysis,  the 
formulation  of  a  model  for  the  temporal  correlation  was 
based  on  past  models  and  ease  of  implementation.  Depending 
upon  the  model  chosen,  the  representation  of  the  temporal 
correlations  required  one  or  more  states  for  each  measure¬ 
ment  pixel,  a  minimum  of  64  new  states  to  add  to  the  truth 
model.  A  first  order  model  was  chosen  because  analysis  had 
not  indicated  the  need  for  a  higher  order  model  and  the  num¬ 
ber  of  states  would  be  minimized.  The  exponential  time  cor¬ 
relation  form  was  assumed  to  be  an  appropriate  model  which 
is  easily  handled  conceptually  and  often  used  in  temporal 
modelling. 

The  stationary  exponentially  correlated  model  of  the 
relationship  between  the  background  noises  of  a  particular 
pixel  at  times  t  and  s  incorporates  the  following  statistics: 

E{v(t)}  =  E(vCs)}  =  '0  (83) 

E{v(t)v(s) }  =  c*  e~'t‘s'/TN  (84) 

where  a„  =  the  desired  variance  of  the  noise 

N 

t„  ■  the  correlation  time  constant  of  the  noise 
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The  continuous  time  equation  for  generating  v(t)  for  a  sin¬ 
gle  pixel  is  represented  in  Figure  13  and  written  as 


v(t)  =  -  1/tn  v  (t )  +  wN(t) 


(85) 


where  w„  is  zero  mean  white  Gaussian  noise  whose  variance 
N 

kernel  must  be  (Ref  8:178) 

2on 

E{wN(t)wN(t+T)}  =  6  (T)  (86) 

N 


so  that  the  variance  of  v(t) 


is  a 


2 

N* 


Figure  13.  Continuous  First  Order  Noise  Generation  Model 

The  equivalent  discrete  time  form  of  this  model  (con¬ 
tinuous  to  discrete  time  conversions  is  covered  in  Ref  8, 
pages  170  to  174)  is  shown  in  Figure  14  and  written  as 


T(ti+i)  -  +  ‘WV 


(67) 


where 


iu^i+i'V  + 


s-(ti+l-ti.)/TN  (88) 


and 


E{wNd(ti)}  "  0 

E{wNd(ti)}  "  °N(1  “  e”2 (ti+l_ti)^TN) 


(89) 

(90) 
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in  order  to  generate  v(ti>  with  appropriate  statistics. 

For  computer  simulation,  the  process  to  generate  v(ti+1)  i 
shown  in  Figure  15.  w(t^)  is  a  zero  mean,  unit  variance, 

white  discrete-time  Gaussian  noise  generated  by  a  Gaussian 
noise  generator. 


wMJ(t.. )  + 


v(t,  . ,  ) 


Delay 


Figure  14.  Discrete  First  Order  Noise  Generation  Model 


»“  (ti+l^i)  /tn1<  — - 1  Delay 


Figure  15.  Discrete  First  Order  Noise  Generation  Model 
Driven  by  Unit  Variance  Noise 


Since  the  64  noise  states  are  independent  of  each  other 
the  equations  for  v(t^+^} ,  a  vector  of  dimension  64,  is 
easily  written  as 

2(ti+i>  -  +  S*s(ti>  (9l> 

where  l^^i+l^i^  *s  a  64-by-64  diagonal  matrix  written  as: 


^N(ti+l'ti) 


s-(ti+i-ti)/TN  j 


The  covariance  of  w^tt^)  is  diagonal  matrix  of  dimensions 
64  by  64 

E{2Nd(ti)sL(<W>  -B-  °N(1  -  e'2(ti+1'ti)/TN)I 

(93) 

Since  these  states  are  decoupled  from  the  other  truth  model 
states  described  in  (54) ,  the  64  noise  states  can  be  aug¬ 
mented  as  a  block  onto  the  end  of  the  original  eight  state 
truth  model,  yielding  a  72  state  truth  model. 

The  generation  of  background  noise  which  is  spatially 
correlated  as  well  as  temporally  correlated  is  accomplished 
by  allowing  the  driving  noise  of  the  64  shaping  filters  to 
be  correlated.  This  is  simply  accomplished  by  inserting  the 
appropriate  nonzero  values  in  the  off-diagonal  terms  in  (93) 
above.  The  modelling  and  computation  of  R  for  spatially 
correlated  noises  is  detailed  in  a  preceding  section  titled 
Spatially  Correlated  Background  Noise.  This  formulation 
assumes  the  separability  of  the  spatial  and  temporal  effects 
in  the  background  noise. 

The  implementation  of  the  temporally  correlated  back¬ 
ground  noises  requires  appropriate  state  rearrangement 
(index  changes)  to  compensate  for  the  apparent  motion  of  the 
background  due  to  target  motion.  When  tracking  a  target 
moving  horizontally  from  left  to  right,  the  background  will 
appear  in  the  FLIR  to  move  with  equal  velocity  in  the  oppo¬ 
site  horizontal  direction  (right  to  left) .  To  keep  the 
simulated  background  of  the  target  correct,  three  things  must 
be  done:  the  target  motion  for  each  sample  period  must  be 
computed,  the  temporal  noise  states  must  be  moved  to  reflect 
the  target  motion,  and  those  states  representing  pixels  with 
new  background  roust  be  filled  with  new  spatially  correlated 
noises  (not  time  correlated  with  previous  state  values 
because  the  previous  history  of  these  pixels  is  not  known) . 
These  three  steps  make  up  three  subroutines  in  the  simula¬ 
tion  software  which  can  be  seen  in  Appendix  F.  The  pixel 
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motion  taken  precisely  from  the  other  truth  model  states  is 
rounded  to  the  nearest  integer  pixel  before  the  state 
rearrangements  are  performed.  The  motion  was  restricted  to 
integer  pixel  values  to  allow  the  same  8-by-8  array  model 
to  be  used  at  all  times  easily.  The  rounded  motion  value 
is  a  fraction  of  a  pixel  in  error  from  the  true  motion 
value.  This  fraction  is  stored  for  one  sample  period  and 
added  to  the  next  motion  value  before  the  rounding  process 
is  executed.  ■  This  avoids  an  accumulation  of  error  due  to 
the  rounding  process.  This  approximation  was  deemed  valid 
for  target  velocities  of  significant  size,  thirty  or  more 
pixels  per  second  (or  for  motionless  targets,  0  pixels  per 
second,  when  no  state  transitions  are  needed) . 

In  summary,  a  step  by  step  algorithm  for  the  generation 
of  temporally  and  spatially  correlated  noises  is  as  follows. 

(1)  w,  a  64-dimensional  vector  of  zero  mean,  unit  var¬ 
iance,  white  Gaussian  noises  is  generated. 

(2)  wNd  is  computed  by  premultiplying  w  with  ^/R r  the 
Cholesky  square  root  of  R. 

(3)  The  states  v  are  propagated  forward  in  time  by 
multiplying  by  ^(t^^t^)  and  then  adding  w^. 

(4)  The  pixel  motion  for  the  time  step  is  computed  to 
the  nearest  pixel  and  the  ignored  fractions  are 
stored . 

(5)  The  states  are  moved  appropriately. 

C6)  Those  states  now  representing  new  background  are 
filled  with  new  noises. 

Field  Size 

The  FLIR  field  of  view  (FOV)  is  a  fixed  angle  and  in 
this  research  was  taken  to  be  160  prads  wide.  Although 
arbitrarily  chosen,  this  is  close  to  the  values  AFWL  is 
using.  As  a  vehicle  moves  toward  or  away  from  the  FLIR,  the 
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angle  subtended  by  the  image  increases  or  decreases.  At 
some  range  the  image  will  fill  the  FOV.  For  example,  at 
10  km  a  vehicle  1.6  m  long  will  fill  the  FOV.  Clearly, 
this  is  inadequate  for  tracking  air-to-air  missiles.  Two 
alternatives  seem  readily  apparent,  a  larger  FOV  or  compres¬ 
sing  the  image  to  stay  within  the  existing  FOV.  That  is, 
combine  a  2  x  2  array  (4  elements)  into  one  pixel,  etc. 

Larger  FOV  seems  attractive  except  that  it  would 
increase  the  computational  burden  significantly.  Increasing 
the  FOV  from  the  current  8x8  pixel  array  to,  say,  16  x  16 
will  quadruple  the  computations  required  in  some  parts  of 
the  filter  algorithm.  When  computation  time  becomes  a 
large  part  of  the  sample  period,  the  assumption  that  control 
outputs  occur  almost  instantaneously  with  sample  time  is  no 
longer  valid.  This  may  introduce  instability  into  a  system 
which  otherwise  would  be  stable. 

Compressing  the  image  will  keep  it  within  the  FOV  but 
at  the  expense  of  resolution.  However,  if  the  tracking  sys¬ 
tem  works  well,  the  loss  in  accuracy  may  be  tolerable.  Two 
compression  simulation  approaches  were  considered:  (1)  have 
the  truth  model  actually  compute  a  16  x  16  array  to  present 
to  the  filter  for  2:1  compression,  or  (.2)  reduce  the  truth 
model  sigma's  by  some  compression  factor  but  still  present 
the  filter  an  8  x  8  array.  Either  way  would  have  required 
extensive  programming  while  not  furthering  the  main  thrust 
of  this  research.  Therefore,  image  dimensions  were  main¬ 
tained  within  the  FOV  to  avoid  the  need  to  confront  the 
compression  issue. 

Although  not  done  here,  consideration  of  FOV  compression 
or  expansion  is  an  area  which  is  a  logical  extension  of  this 
research . 


Referencing 


REFERENCING,  which  is  a  tracking  simulation  package 
supplied  by  AFWL,  was  used  primarily  for  analyzing  the 
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temporal  correlations  of  the  background  noise  in  real  data 
as  described  previously.  However,  REFERENCING  has  consid¬ 
erable  potential  to  be  used  as  an  alternate  truth  model  and 
as  a  more  extensive  analysis  tool. 

REFERENCING  has  two  operating  modes:  one  mode  uses 
real  digital  target  data  from  magnetic  tape  and  the  other 
mode  generates  its  own  synthetic  target  data.  In  either 
mode  the  real  target  or  synthetic  target  is  tracked  by  one 
of  six  correlation  tracking  schemes  currently  being  studied 
by  AFWL  CRef  13:1-2).  Starting  with  a  16-by-16  reference 
array  containing  an  edge  of  the  target  from  an  initial  data 
frame,  REFERENCING  computes  a  16-by-16  (or  smaller)  matrix 
of  cross-correlations  between  the  reference  array  and  sub¬ 
matrices  of  a  search  matrix  in  a  subsequent  data  frame. 

(The  cross  correlation  array  was  used  in  temporal  noise 
analysis  in  Appendix  C.)  The  cross  correlation  matrix  is 
then  utilized  by  the  correlation  trackers  to  compute  a  pre¬ 
cise  target  offset  between  the  two  frames. 

Compared  to  the  truth  model  generated  in  this  study, 
REFERENCING  has  one  serious  drawback  as  a  good  representa¬ 
tion  of  the  real  world.  In  terms  of  target  size  and  shape, 
and  the  noise  characteristics  of  the  measurements,  the  real 
digital  data  from  tape  cannot  be  surpassed  by  a  software 
model  of  the  target  and  noises.  However,  because  the  real 
data  was  taken  from  the  FLIR  which  was  moving  to  keep  the 
target  centered  in  its  field  of  view,  the  target  appears 
essentially  motionless.  Moreover,  the  true  motion  (states) 
of  the  target  is  not  discernible  from  the  data  tapes.  Thus, 
the  errors  committed  by  any  tracking  system  could  not  be 
computed . 

When  operated  in  the  synthetic  target  mode,  REFERENCING 
is  deficient  in  representing  any  real  target  motion.  The 
synthetic  target  can  be  generated  in  almost  any  size  and 
shape  with  its  edges  smoothed  (Gaussian  shape)  by  a  filter; 
in  terms  of  target  intensity  characteristics  it  is  a  good 


simulation.  The  target  can  then  be  either  segmented  (rows 
or  columns  of  the  target  deleted) ,  rotated,  or  jittered 
during  each  run.  However,  the  target  cannot  display  real 
motion.  Thus,  the  target  motion  is  again  lacking  in  real¬ 
ism,  and  except  when  the  target  is  jittered  there  is  no 
noise  in  the  synthetic  target  frames. 

With  help  from  AFWL,  the  REFERENCING  routine  could  be 
made  into  an  excellent  "truth  model"  using  real  digital 
data.  The  target  motion  could  be  recreated  from  the  time 
history  of  the  FLIR  pointing  angles  and  any  other  target 
data  available  (i.e.  radar  range  and  velocity)  from  AFWL. 
The  true  trajectory  could  not  be  computed,  but  an  adequate 
noise-corrupted  representation  probably  could  be  achieved. 


III.  Previous  Filter  Robustness 


Introduction 

The  robustness  of  the  original  filter  is  investigated 
when  it  is  subjected  to  truth  model  parameter  mismatches 
and  to  major  truth  model  changes.  The  first  group  of  com¬ 
puter  simulations,  using  the  original  truth  model  and  fil¬ 
ter,  investigated  the  filter  performance  when  parameter 
mismatches  existed  between  the  truth  model  and  the  filter. 

A  second  group  of  simulations  evaluated  the  performance  of 
the  original  filter  when  subjected  to  the  new  target  sce¬ 
nario  and  background  as  depicted  by  major  truth  model 
changes  described  in  Chapter  II. 

Before  discussing  the  particular  computer  simulation 
runs,  the  general  simulation  scenario  will  be  briefly  out¬ 
lined.  Monte  Carlo  simulations  were  performed  on  a  CDC  6600 
using  software  packages  written  in  FORTRAN  IV.  Each  Monte 
Carlo  case  consisted  of  20  five  second  simulations.  Each 
of  the  20  simulations  was  started  at  the  same  specific 

value  of  x  (0)  and  P  (.51).  The  data  from  the  Monte  Carlo 
— o  —  o  — 

runs  was  later  processed  and  appropriate  statistics  for  the 
atmospheric  jitter  and  target  state  estimate  errors  were 
plotted  on  graphs.  Appendix  E  discusses  in  more  detail 
notable  computer  simulation  considerations  such  as  noise 
generation,  approximations,  and  tuning.  Appendix  F  con¬ 
tains  copies  of  the  FORTRAN  IV  software  packages  used  to 
perform  the  Monte  Carlo  simulation.  Appendix  G  contains  a 
copy  of  the  software  module  which  computed  the  simulation 
statistics  and  plotted  the  data.  Appendix  H  lists  the 
parameters  for  cases  1-31  which  were  performed  as  part  of 
the  analysis  for  this  chapter.  Appendices  I  and  J  contain 
the  plots  which  display  the  detailed  results  of  the  simula¬ 
tions  contained  in  this  chapter;  major  results,  trends,  and 
discussions  will  be  presented  in  the  chapter  itself. 
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Sensitivity  to  Parameter  Mismatches 


The  performance  of  the  previously  designed  filter  was 
evaluated  by  inputting  contrasting  parameters  for  the  filter 
and  truth  models.  Two  preliminary  baseline  cases  were  run. 
Case  1  was  run  at  a  signal  to  noise  ratio  of  20  (S/N=I  Jo„) 
and  served  as  a  comparison  to  previously  reported  results 
(Ref  5) .  Case  2  was  run  as  a  baseline  for  cases  where 
S/N=10.  With  a  nominal  set  of  parameters,  Case  2  serves  a 
performance  standard  for  the  other  runs.  Cases  3  through  12 
duplicate  case  2  except  one  parameter  in  the  truth  model  was 
changed  in  each  case.  Five  areas  of  sensitivity  were  exam¬ 
ined  by  giving  the  truth  model  high  and  low  values  of  a 
parameter  while  the  filter  was  given  the  nominal  value  of 
the  parameter.  The  areas  of  interest  were  target  dynamics 
correlation  time,  target  intensity,  target  size,  the  root 
mean  square  (rms)  value  of  background  noise,  and  the  ratio 
of  the  target  dynamics  rms  value  to  the  atmospheric  jitter 
rms  value. 

Case  1,  whose  input  parameters  can  be  found  in  Appendix 
H  and  resulting  plots  can  be  found  in  Appendix  I,  was  run 
to  check  out  the  inherited  software  packages  through  com¬ 
parison  to  Case  1  conducted  by  Captain  Mercier  (Ref  5:47). 
Both  cases  were  run  at  a  signal  to  noise  (S/M)  ratio  of  20, 
and  they  agreed  very  closely  in  final  results  (see  Appendix 
I  and  Ref  5:96  to  98).  Mercier' s  Case  1  was  not  duplicated 
exactly  due  probably  to  a  difference  in  the  seeds  for  the 
random  number  generator.  However,  the  variances  of  the  real 
error  and  of  the  filter  estimated  error  agreed  very  closely. 

Case  2,  which  was  conducted  at  a  signal  to  noise  ratio 
of  10,  serves  as  a  baseline  for  comparison  for  the  rest  of 
the  mismatch  cases.  Case  2  was  executed  with  all  the  nomi¬ 
nal  parameters,  no  mismatched  parameters.  Table  I  summa¬ 
rizes  the  results  of  case  2  -  the  units  are  pixels  for  the 
average  mean  absolute  error  and  for  the  average  standard 
deviation.  Because  all  the  errors  tend  to  be  zero  mean, 


TABLE  I 


Baseline  for  Cases  3-18: 

Case  2 

Results 

X 

y 

Average  Standard  Deviation 

Atmospheric  Jitter 

Actual 

.55 

.58 

Filter-Indicated 

.58 

.58 

Target  Position 

Actual 

.56 

.58 

Filter-Indicated 

.58 

.58 

Average  Absolute  Error 

Atmospheric  Jitter 

.07 

.08 

Target  Position 

.10 

.12 

the  average  mean  absolute  error  proved  to  be  more  meaning¬ 
ful  and  is  presented  in  the  tables  to  follow.  Each  value 
shown  represents  the  average  of  the  value  over  the  five 
second  simulation  period.  Figures  16,  17,  and  18  are  three 
of  the  plots  resulting  from  case  2.  The  three  plots  display 
statistics  for  the  target  dynamics  in  the  x  direction. 
Similar  plots  were  made  for  the  x-direction  atmospheric  jit¬ 
ter,  the  y-direction  target  dynamics,  and  the  y-direction 
atmospheric  jitter;  they  can  be  seen  in  Appendix  I.  Figure 
16  demonstrates  the  true  error  standard  deviation  (0)  con¬ 
vergence  of  the  target  dynamics  estimate.  The  true  error 
standard  deviation  (0)  is  plotted  as  a  function  of  the  num¬ 
ber  of  runs  to  demonstrate  the  convergence  which  occurs 
when  20  simulations  are  run.  Figure  17  is  a  very  effective 
tool  for  tuning  analysis  as  it  displays  the  true  error  stan¬ 
dard  deviation  (the  noisy  curve)  as  well  as  the  filter- 
indicated  error  standard  deviation  (the  more  consistent 
curve) ,  good  tuning  being  achieved  typically  when  these 
agree.  Figure  18  displays  the  mean  true  target  position 
error  for  the  20  Monte  Carlo  runs  and  the  one  sigma  envelope 
of  the  errors  for  performance  analysis. 
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Figure  16.  Case  2:  Target  Position  x  Channel  Standard 
Deviation  (a,  Sigma)  Convergence 
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Figure  17.  Case  2:  Target  Position ,  X  Channel,  Real  Error 
versus  Filter  Standard  Deviation 
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The  first  area  of  interest  for  studying  filter  robust¬ 
ness  was  variations  in  target  dynamics.  The  nominal  value 
for  td  in  equations  C3)  and  (4)  was  set  at  1  sec.  For  case 
3  the  filter  assumed  td  =  1,  but  the  truth  model  was  based 
on  Tjj  =  0.2.  This  can  be  interpreted  as  the  target  exhibit¬ 
ing  higher  frequency  motion  (more  jittery)  than  the  filter 
assumes.  For  case  4,  the  truth  model  is  based  on  =  5,  a 
slower,  more  time-correlated  motion  than  the  filter  expects. 
The  results  for  both  cases  are  shown  in  Table  II  and  a  lim¬ 
ited  set  of  their  plots  can  be  seen  in  Appendix  I.  The  fil¬ 
ter  performance  for  both  cases  is  very  acceptable.  However, 
the  results  have  little  meaning  in  terms  of  the  air-to-air 
missile  scenario  because  the  truth  model's  target  motion 
equation  for  these  cases  is  extremely  benign. 

Cases  5  and  6  test  the  filter's  response  to  mismatches 
in  maximum  target  intensity.  The  nominal  maximum  target 
intensity  (ITn-,v)  measure  of  10  (intensity  units,  the  exact 
units  are  irrelevant  to  the  purpose  of  this  as  long  as 
appropriate  ratios  are  maintained)  was  used  by  the  filter 
with  the  rms  value  of  the  background  noise  set  at  1  (intens¬ 
ity  unit).  For  case  5  the  truth  model  is  given  I  =  1; 
the  real  FLIR  images  show  a  target  image  that  is  very  much 
masked  by  the  noise  while  the  filter  thinks  that  the  target 
image  should  be  readily  distinguishable.  The  truth  model 

for  case  6  has  I  =20.  Now  the  target  intensity  is 
max 

greater  than  the  filter  model  assumes.  The  performance  of 
the  filter  in  both  cases  is  not  acceptable  as  can  be  seen 
in  Table  III  and  their  plots  in  Appendix  I.  While  the  poor 
performance  of  case  5  is  hardly  surprising  (such  a  low  S/N 
can  cause  poor  performance  with  or  without  a  mismatch) ,  the 
fact  that  case  6  exhibited  similar  results  is  very  interest¬ 
ing.  The  filter  appears  to  have  difficulty  because  it  is 
searching  for  the  wrong  shape  of  intensity  profile  caused  by 

the  mismodelled  I  .  These  results  should  generally  be 

max 

applicable  to  the  close  range  air-to-air  missile  problem 
because  both  scenarios  share  the  Gaussian  intensity  model. 


TABLE  II 


Target  Dynamics  Mismatch  - 
Case  3:  Low  xDT  and  Case  4:  High  tDT  Results 


Case 

3 

Case 

4 

TDT  = 

.2 

TDT  ■ 

5 

Average  Standard  Deviation 

X 

y 

X 

y 

Atmospheric  Jitter 

Actual 

.65 

.65 

.48 

.45 

Filter-Indicated 

.58 

.58 

.58 

.58 

Target  Position 

Actual 

.65 

.68 

.45 

.42 

Filter-Indicated 

.58 

.58 

.58 

.58 

Average  Mean  Absolute  Error 

Atmospheric  Jitter 

.10 

.12 

.10 

.07 

Target  Position 

.11 

.11 

.11 

.09 

The  third  set  of  mismatches  that  were  studied  dealt 
with  target  size.  The  nominal,  filter,  value  for  the  target 
size  was  three  pixels  (ag  =  3) .  In  case  7  the  truth  model 
had  the  target  size  set  at  1  pixel  (a  =  1) ;  the  real  target 
is  much  smaller  than  the  target  the  filter  is  trying  to 
identify.  Consequently,  the  real  image  randomly  walks  in 
the  large  target  envelope  that  the  filter  is  searching  for, 
with  the  result  of  poor  performance,  as  can  be  seen  in 
Table  IV  and  the  plots  in  Appendix  I.  The  filter  performs 
acceptably  in  case  8,  in  which  **  5  for  the  truth  model. 
The  filter  appears  to  pin  down  the  center  of  the  target  much 
better  if  the  real  image  is  larger  than  the  assumed  filter 
image  size,  as  can  be  seen  in  Table  IV.  These  results  are 
very  meaningful  for  the  close  range  target  where  the  real 
image  size  may  be  not  only  mismatched  but  also  changing  in 
size  and  shape. 

Mismatches  in  the  background  noise  were  studied  in 
cases  9  and  10.  The  nominal  value  for  the  signal  to  noise 
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TABLE  III 


Target  Intensity  Mismatch  - 
Case  5:  Low  Imax  and  Case  6:  High  Imax  Results 


Case  5 

Case 

6 

I 

max 

=  1 

Imax 

>  20 

Average  Standard  Deviation 

X 

y 

x 

y 

Atmospheric  Jitter 

Actual 

1.7 

1.60 

1.5 

1.5 

Filter-Indicated 

.60 

.60 

.58 

.58 

Target  Position 

Actual 

3.5 

4.0 

.70 

.70 

Filter-Indicated 

.65 

.65 

.58 

.58 

Average  Absolute  Error 

Atmospheric  Jitter 

.20 

.40 

.25 

.25 

Target  Position 

.20 

.10 

.10 

.20 

ratio  (S/N  =  Imav/gH)  was  10,  as  modelled  in  the  filter. 

For  case  9  the  truth  model  S/N  was  1  (the  variance  of  the 
background  noise  is  10  instead  of  1  and  IJnax  equals  10  for 
both  filter  and  truth  model) .  This  case  is  very  similar 
to  case  5;  the  filter  is  searching  for  a  prominent  target 
against  low  background  noise  but  the  real  target  is  actually 
surrounded  by  background  noise  of  equal  intensity.  In  case 
10,  the  truth  model  has  a  S/N  =  20  so  that  the  variance  of 
the  background  noise  is  0.5.  Case  9  resulted  in  terrible 
filter  performance  (Table  V)  for  reasons  similar  to  those 
of  the  poor  performance  exhibited  in  case  5.  For  both 
cases,  the  real  measurements  in  both  portray  a  target  image 
effectively  hidden  among  the  noise.  This  is  probably  the 
major  cause  of  the  poor  performance  in  both  cases  rather 
than  the  mismatches  which  existed.  A  similar  analogy  does 
not  exist  between  cases  6  and  10.  Case  10  exhibited  very 
good  filter  performance  which  is  understandable  when  inter¬ 
preted  physically.  Both  the  truth  model  and  filter  in  case 

10  have  I  „„  ®  10.  The  mismatch  is  that  the  S/N=10  for  the 
max 


53 


TABLE  IV 


Target  Size  Mismatches  - 


Case  7:  a  =  1  and  Case 
gT 

8: 

°gT  ~ 

5  Results 

Case 

7 

Case 

8 

• 

°gT 

=  1 

°gT  = 

5 

Average  Standard  Deviation 

X 

y 

X 

y 

Atmospheric  Jitter 

Actual  1 

.5 

1.5 

.63 

.62 

Filter-Indicated 

.58 

.58 

.58 

.58 

Target  Position 

Actual  3 

.0 

3.0 

.65 

.65 

Filter-Indicated 

.61 

.62 

.58 

.58 

Average  Absolute  Error 

Atmospheric  Jitter 

.25 

.4 

.10 

.12 

Target  Position 

.5 

1.5 

.15 

.10 

filter  results  in  the  variance  of 

the 

noise 

being  1 

whereas 

the  variance  of  the  noise  in  the  truth  model  is  0.5  (S/N=20). 
Simply,  the  assumed  target  intensity  shape  is  correct 
(instead  of  erroneous  as  in  case  6)  while  the  background 
noise  is  actually  less  than  the  filter  assumes,  and  it  tracks 
the  target  very  well. 

The  final  mismatch  runs  investigated  the  ratio  of  rms 
value  of  the  target  motion  to  rms  value  of  the  atmospheric 
jitter.  These  values  are  assumed  to  be  equal  for  the  nominal 
case  (oD/oA  =1).  In  case  11,  the  rms  value  of  the  target 
motion  was  set  to  a  fifth  of  the  rms  value  of  atmospheric 
jitter  (^D/^A  =  .2).  Physically  interpreted,  the  real  tar¬ 
get  motion  is  much  less  in  magnitude  than  the  filter  is 
searching  for.  So  the  filter  over  estimates  the  errors, 
as  can  be  seen  in  Table  VI.  In  case  14  the  real  target  has 
considerably  more  motion  taD/tfA  =  5}  than  the  filter  is  pro¬ 
grammed  to  expect,  resulting  in  unacceptable  filter  perform¬ 
ance  shown  again  in  Table  VI  and  Appendix  I.  The  results 
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TABLE  V 


Background  Noise  Mismatch  - 


Case  9:  S/N  =  1  and  Case 

10: 

S/N  = 

20  Results 

Case  9 

Case  10 

• 

S/N 

=  1 

S/N 

=  20 

Average  Standard  Deviation 
Atmospheric  Jitter 

X 

y 

X 

y 

Actual 

8.0 

8.1 

.57 

.56 

Filter-Indicated 

Target  Position 

.63 

.63 

.58 

.58 

Actual 

3.8 

3.9 

.56 

.53 

Filter-Indicated 

Average  Absolute  Error 

.77 

.78 

.58 

.58 

Atmospheric  Jitter 

1.5 

1.5 

.12 

.10 

Target  Position 

.5 

.6 

.12 

.12 

these  cases  again  have  little 

meaning  in 

terms 

of  the 

air-to-air  missile  problem.  The  truth  model's  target  motion 
and  the  filter's  target  motion  model  in  these  cases  are 
totally  unrealistic  in  terms  of  representing  an  air-to-air 
missile.  The  truth  model's  target  motion  has  already  been 
redefined  in  Chapter  II.  An  appropriate  change  in  the 
filter's  target  dynamics  model  is  therefore  motivated. 

Sensitivity  to  New  Truth  Model 

The  previously  formulated  filter  was  also  evaluated  as 
it  was  subjected  to  the  major  truth  model  changes  described 
in  Chapter  II.  First,  the  filter  was  tested  against  targets 
with  differing  aspect  ratios  CAR  =  image  length/image  width 
**  <Jv/OpV  as  defined  in  Chapter  II)  as  simulated  by  the  truth 
model.  Next,  filter  performance  was  evaluated  when  the 
background  noise  was  spatially  correlated.  The  truth  model 
was  then  changed  to  allow  target  motion,  and  the  filter 
performance  tested  in  terms  of  the  new  motion.  Temporally 
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TABLE  VI 


Ratio  of  Atmospheric  Jitter  to  Target  Dynamics  Mismatch  - 
Case  11:  oD/oA  =  *2  and  Case  12:  °D/°A  =  5  Results 


Case 

11 

Case 

12 

Average  Standard  Deviation 
Atmospheric  Jitter 

X 

y 

X 

y 

Actual 

.43 

.41 

1.7 

1.5 

Filter-Indicated 

Target  Position 

.58 

.58 

.58 

.58 

Actual 

.38 

.37 

2.5 

2.8 

Filter-Indicated 

.58 

.58 

.59 

.60 

Average  Absolute  Error 

Atmospheric  Jitter 

.08 

.08 

.35 

.4 

Target  Position 

.10 

.07 

.4 

.5 

correlated  background  noise  was  the  last  addition  to  the 
truth  model,  and,  subsequently,  the  filter  performance  was 
evaluated  for  various  target  scenarios  with  a  variety  of 
temporal  and  spatial  background  noises.  Each  truth  model 
change  as  formulated  in  Chapter  II  was  aimed  at  a  more 
realistic  representation  of  the  problem.  The  filter  per¬ 
formance  in  light  of  these  changes  to  the  truth  model  was 
not  expected  to  be  good  because  the  filter  had  not  been 
designed  to  track  the  type  of  target  which  the  changes  were 
attempting  to  create.  The  goal  of  these  runs  was  to  obtain 
a  prioritized  list  of  those  characteristics  of  the  new  tar¬ 
get  scenario  which  caused  the  most  problems. 

The  filter's  ability  to  track  a  target  with  an  aspect 
ratio  CAR)  greater  than  one  was  first  tested.  This  was 
accomplished  by  a  software  change  to  replace  the  circular 
constant  intensity  contour  target  by  a  target  (also  bivari¬ 
ate  Gaussian)  with  variable  minor  to.)  and  major  Co...  *AR) 

pv  pv 

axes'  magnitudes;  target  dynamics  and  atmospheric  jitter 
state  equations  in  the  truth  model  remained  as  originally 
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formulated.  For  the  cases  tested,  the  major  axis  was 
aligned  with  the  x-axis.  Cases  13,  14,  and  15  were  run 
with  0pv  equal  to  1  and  AR  equal  to  2,  5,  and  10,  respec¬ 
tively.  The  results  of  these  cases  are  shown  in  Table  VII. 
(The  plots  and  input  parameters  for  the  truth  model  and 
filter  for  each  case  are  found  in  Appendices  H  and  I, 
respectively.)  As  expected  in  all  cases,  the  y-axis  filter 
performance  was  very  good  because  the  target  size  did  not 
change  in  that  direction.  In  the  x-axis  direction,  the 
filter  performance  did  not  degrade  significantly  as  the 
target  increased  in  aspect  ratio.  These  results  show  small 
but  definite  errors  which  would  be  amplified  by  other  com¬ 
plicating  factors  such  as  feedback,  target  rotation,  target 
motion,  and  variations  in  target  size.  Other  factors  and 
simulations  need  to  be  considered  before  judging  the  role  of 
AR. 

The  truth  model  changes,  discussed  in  Chapter  II,  for 
creating  spatially  correlated  background  noise  were  imple¬ 
mented,  and  the  filter  was  tested  against  it.  Again,  the 
state  equations  were  not  changed,  and  the  circular  constant 
intensity  contour  target  was  used  (AR  =  1) .  First  and  sec¬ 
ond  spatial  correlations  using  the  exponential  model  were 
used.  The  24  correlation  coefficients  of  any  one  pixel 
have  only  five  different  values  due  to  the  symmetry  of  the 
correlations.  For  case  16  for  example,  r^  ^+1  rjc  ^+9* 

rk ,k+2 *  rk,k+10 '  and  rk,k+18  shown  in  Pi^e  5  were  set  to 
.368,  .243,  .135,  .105,  and  .059,  respectively,  as  generated 
by  an  exponential  correlation  model.  The  other  19  coeffi¬ 
cients  were  set  equal  to  one  of  these  values  as  dictated  by 
the  symmetry.  Cases  17  and  18  were  run  with  spatial  corre¬ 
lations  of  .497,  .372,  .247,  .209,  and  .138  at  S/N  and  Imav 
equal  to  10  and  1,  respectively.  These  coefficients  repre¬ 
sent  the  upper  limit  of  values  which  could  be  used  without 
numerical  problems  in  R  (see  Chapter  II) . 
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TABLE  VII 


Elliptical  Target  Intensity  Contours  in  the  Truth  Model 
Cases  13/  14,  and  15  Results 


Case 

13 

Case 

14 

Case 

15 

AR  = 

2 

AR  * 

5 

AR  = 

10 

X 

y 

X 

y 

X 

y 

Average  Standard  Deviation 

Atmospheric  Jitter 

Actual 

.65 

.58 

.82 

.58 

.85 

.56 

Filter- Ind i ca ted 

.58 

.58 

.58 

.58 

.58 

.58 

Target  Position 

Actual 

.60 

.55 

.80 

.58 

.90 

.54 

Filter- Indicated 

.58 

.58 

.58 

.58 

.58 

.58 

Average  Absolute  Error 

Atmospheric  Jitter 

.15 

.10 

.20 

.12 

.20 

.12 

Target  Position 

.10 

.10 

.15 

.10 

.15 

.10 

Note:  the  baseline, 
is  summarized 

Case  2,  for  comparison 
in  Table  I. 

Cases  16 ,  17 ,  and  18  can  be  interpreted  physically  as 
tracking  a  nearly  stationary  target  against  a  spatially 
correlated  background  while  the  filter  assumes  that  the 
background  noise  is  spatially  white.  The  results  of  these 
three  cases,  which  are  displayed  in  Table  VIII,  showed 
acceptable  performance  except  for  case  18  when  S/N  =  1.  The 
filter  performance  in  cases  16  and  17  was  almost  identical 
to  case  2;  the  spatially  correlated  background  noise  caused 
no  perceived  degradation.  Case  18' s  poor  performance  was 
probably  caused  equally  by  the  low  signal  to  noise  ratio 
and  the  spatial  correlations  in  the  background  noise.  This 
was  partially  confirmed  by  comparison  to  case  10  of  Capt. 
Mercier’s  thesis  (Ref  5:124-125)  which  shows  a  true  error 
variance  for  the  atmospheric  jitter  and  target  position  of 
1.0  and  1.25,  respectively  (even  here  the  filter  underesti¬ 
mated  its  own  errors) .  This  case  corresponds  to  case  18 


TABLE  VIII 


Spatially  Correlated  Background  Noise  in  the  Truth  Model 
Cases  16,  17,  and  18  Results 


Case  16  Case  17  Case  18 


. 

S/N 

=  10 

S/N  =  10 

S/N  =  1 

Low 

Cor- 

High  Cor¬ 

High  Cor¬ 

relations 

relations 

relations 

X 

y 

x  y 

x  y 

Average  Standard  Deviation 

Atmospheric  Jitter 

Actual 

.58 

.57 

.58  .58 

1.5  1.55 

Filter- Indicated 

.58 

.58 

.58  .58 

.85  .85 

Target  Position 

Actual 

.55 

.56 

.58  .56 

2.5  2.2 

Filter- Indicated 

.58 

.58 

.58  .58 

.65  .65 

Average  Absolute  Error 

Atmospheric  Jitter 

.13 

.12 

.12  .10 

.30  .30 

Target  Position 

.08 

.10 

H 

H 

• 

00 

o 

• 

.40  .50 

Notes  the  baseline, 
is  summarized 

Case  2 ,  for 
in  Table  I. 

comparison 

here,  but  with  no  spatial  and  temporal  correlation  in  the 
truth  model  background  noise,  i.e.  with  no  filter  mismatch. 

The  next  test  for  the  filter  was  the  addition  of  motion 
to  the  truth  model  (refer  to  Chapter  II  for  description  of 
target  motion  formulation) .  The  elliptical  target  and  spa¬ 
tial  background  noise  were  not  used  for  these  test  cases. 

One  change  was  made  to  the  filter  to  allow  it  to  try  to 
track  a  moving  target.  The  original  filter  was  open  loop 
in  design;  it  was  assumed  that  the  computed  offsets  were  not 
zeroed  out  by  the  system's  control  mechanism.  Thus,  the 
filter  relied  on  the  zero  mean  characteristic  of  the  origi¬ 
nal  truth  model's  dynamics  to  keep  the  target  in  the  FLIR 
field  of  view.  Without  implementing  feedback,  the  filter 
would  lose  any  moving  target  in  one  or  two  sample  periods 


depending  on  the  magnitude  of  the  motion.  (For  more  detailed 
analysis  see  Chapter  IV.)  Thus,  the  filter  was  changed  to 
reflect  control  system  zeroing  of  the  computed  target  motion 
offsets  for  each  sample  period.  The  velocities  for  cases  21 
and  22  were  56.4  pixels  per  second  (x  «*  40  pixels/sec, 
y  =  40  pixels/sec,  initially)  and  84.6  pixels  per  second 
(x  *  60  pixels/sec,  y  «  60  pixels/sec,  initially),  respec¬ 
tively. 

Because  of  the  software  change  to  provide  system  feed¬ 
back,  a  new  baseline  case  was  computed.  The  parameters  and 
plots  for  case  19  can  be  found  in  Appendices  H  and  J,  respec¬ 
tively.  Case  19  was  run  with  a  signal  to  noise  ratio  of  10 
and  without  target  motion  (x  =  0,  y  =  0).  The  results  of 
case  19  are  summarized  in  Table  IX.  The  results  show  the 
filter  to  be  slightly  out  of  tune  (.tuning  is  discussed  more 
completely  in  Appendix  E) .  Tuning  was  not  accomplished 
because  of  time  constraints  and  the  small  amount  the  filter 
was  mis tuned.  This  is  partially  because  of  the  feedback 
change  and  partially  due  to  a  change  in  two  of  the  nominal 
parameters  (oDT/°AT  -  5  and  a ^  =  .2).  Case  19  will  be  the 
baseline  of  comparison  for  cases  20  to  23. 

Table  IX  and  subsequent  tables  in  this  chapter,  a  new 
analysis  value  has  been  added  -  the  number  of  times  in  the 
20  simulations  that  the  target  was  lost.  If  the  true  target 
position  ever  was  more  than  three  times  the  target  size 
C3*av)  away  from  the  filter  computed  position,  the  target 
was  considered  lost.  For  targets  where  cv  is  greater  than 
or  equal  to  1.5  pixel,  this  means  most  of  the  target  image 
has  left  the  field  of  view.  For  any  size  target  this  cri¬ 
teria  means  that  the  true  image  is  so  far  from  the  predicted 
image  that  the  residuals  provide  little  useful  information. 

In  most  cases  where  the  target  is  lost  more  than  one  or  two  - 
times,  the  total  number  of  loss  of  track  is  more  meaningful 
than  the  statistics  that  result. 
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TABLE  IX 


Baseline  for  Cases  20- 

23 

Case  19  Results 

(s/n=i 0) 

Case 

19 

X 

y 

Average  Standard  Deviation 

Atmospheric  Jitter 

Actual 

.25 

.25 

Filter-Indicated 

.19 

.19 

Target  Position 

Actual 

.40 

.40 

Filter-Indicated 

.19 

.19 

Average  Absolute  Error 
Atmospheric  Jitter 

.09 

.08 

Target  Position 

.10 

.10 

Loss  of  track  -  0  times 


The  filter  performance  for  cases  20,  21,  and  22,  sum¬ 
marized  in  Table  X,  show  a  trend  of  increasingly  poor  per¬ 
formance  corresponding  to  the  increasing  velocity  of  the 
cases.  Case  22  statistics  are  not  shown  because  the  loss  of 
track  problems  caused  numerical  problems  serious  enough  to 
interrupt  its  execution.  It  is  known  that  loss  of  track 
occurred  in  all  20  simulations  for  case  22.  The  results  of 
these  three  cases  demonstrate  all  too  well  the  need  for 
better  filter  modelling  of  the  target  motion. 

The  last  truth  model  change  was  the  formulation  of  a 
72  state  truth  model  to  allow  temporal  noise  modelling  of 
the  background  noise.  The  filter  performance  was  analyzed 
in  cases  of  various  combinations  of  S/N,  temporal  correla¬ 
tions,  and  spatial  correlations.  To  make  sense  physically, 
temporal  and  spatial  correlations  were  used  together  in  the 
background  noise  for  most  cases. 

The  first  two  cases,  cases  23  and  24,  were  run  at  S/N 
equal  to  10  and  1,  respectively,  with  both  spatial  and 
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TABLE  X 


Target  Motion  in  the  Truth  Model 
Case  20  and  21  Results 


Case 

20 

Case  21 

x  =  20 

P/S 

x  = 

40  P/S 

y  -  20 

P/S 

y  = 

40  P/S 

X 

y 

X 

y 

Average  Standard  Deviation 

Atmospheric  Jitter 

Actual 

.26 

.28 

.22 

.23 

Filter-Indicated 

.19 

.19 

.19 

.19 

Target  Position 

Actual 

.45 

.50 

8.0 

16.0 

Filter-Indicated 

.20 

.20 

.20 

.20 

Average  Absolute  Error 

Atmospheric  Jitter 

.12 

.14 

.11 

.15 

Target  Position 

.20 

.22 

7.0 

15.0 

Loss  of  Track 

0  times 

12 

times 

Note:  1)  Case  22  lost  track  20  times. 

2)  The  baseline.  Case  19,  for 
comparison  is  summarized  in 
Table  IX. 

temporal  correlations.  For  each  case  the  spatial  correla¬ 
tion  coefficients  were  the  same  as  those  used  in  case  16; 
the  truth  model  depicted  a  motionless  target.  For  the  tem¬ 
poral  noise  tn  was  chosen  (t^  =  20  sec)  to  simulate  a  .95 
correlation  between  noises  over  one  second  of  time.  This 
value  of  tn  is  probably  a  realistic  upper  limit  (worst  case) 
on  the  correlations;  a  correlation  of  1  is  the  theoretical 
upper  limit.  Filter  performance  in  case  23  shown  in  Table 
XI  was  only  slightly  worse  than  that  of  the  baseline,  case 
19,  and  acceptable.  The  performance  of  case  24  was  so  bad 
that  numerical  problems  halted  the  simulation  software. 

These  two  cases  reaffirmed  past  conclusions:  at  S/N  equal 


TABLE  XI 


Temporal  and  Spatial  Background  Noise  in 
Truth  Model  (SN  =  10)  -  Case  23  Results 


Case  23 
•95  Time  Cor. 


Average  Standard  Deviation 
Atmospheric  Jitter 

X 

y 

Actual 

.26 

.27 

Filter-Indicated 

Target  Position 

.20 

.20 

Actual 

.50 

.49 

Filter-Indicated 

.19 

.19 

Average  Absolute  Error 

Atmospheric  Jitter 

.09 

.08 

Target  Position 

.15 

.12 

Loss  of  Track 

0 

times 

Note:  1)  Case  24  loss  of 

track  »  20 

times . 

2)  The  baseline.  Case  20,  for 
comparison  is  summarized  in 
Table  IX. 


to  10,  no  matter  what  realistic  model  used,  the  background 
noise  has  no  significant  effect  on  filter  performance;  at 
S/N  equal  to  1  real  tracking  problems  occur. 

A  group  of  cases,  simulated  with  S/N  equal  to  2,  pro¬ 
vided  more  interesting  evaluations  of  the  temporal  and  spa¬ 
tial  background  noise  correlations  effect.  A  new  baseline, 
case  25,  was  computed  because  of  the  new  S/N.  Case  25  was 
run  without  truth  model  motion  and  with  temporal  and  spatial 
correlations  in  the  background  noise  deactivated  (actually 
tn  •  .01  sec  which  corresponds  to  a  .03  correlation  in  one 
sample  period).  Case  25  results  are  tabulated  in  Table  XII. 
While  the  performance  statistics  are  not  good,  loss  of  track 
only  occurred  once,  and  these  values  will  serve  for  compari¬ 
son  purposes  for  the  cases  to  follow. 
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TABLE  XII 


Baseline  for  Cases  26-29:  Case  25  Results 

(S/N=2) 


Average  Standard  Deviation 

Atmospheric  Jitter 
.  Actual 

Filter-Indicated 

Target  Position 
Actual 

Filter-Indicated 

Average  Absolute  Error 
Atmospheric  Jitter 
Target  Position 


Case  25 


x  y 


.20  .19 

.19  .19 


1.9  1.3 

.34  .34 

.05  .05 

.40  .30 


Loss  of  Track 


1  time 


Cases  26,  27,  28,  and  29  were  run  at  S/N  equal  to  2 
and  various  combinations  of  temporal  and  spatial  background 
noise.  Case  26  had  high  temporal  correlations  (.95  for  one 
second)  and  spatial  correlations  with  the  same  coefficients 
as  cases  17,  24,  and  25.  Case  28  ran  with  no  temporal  cor¬ 
relations  (xN  =  .01)  but  again  with  the  same  spatial  corre¬ 
lations  as  cases  27  and  28 .  Case  29  was  the  reverse  of 
case  28;  no  spatial  correlations  but  high  temporal  correla¬ 
tions  were  used  (.95  for  one  second).  The  results  for  all 
cases  are  shown  in  Table  XIII.  Case  26  portrays  a  target 
against  a  background  with  considerable  spatial  and  time  cor¬ 
related  infrared  structure.  Filter  performance  was  poor 
with  loss  of  track  resulting  in  12  simulations.  Case  27  is 
the  same  as  case  26  but  with  less  time  correlation  in  the 
background  structure.  Performance  was  similarly  poor. 

Case  28  emulates  a  background  that  is  correlated  in  space 
but  not  in  time;  performance,  while  not  acceptable,  was  very 
close  to  that  of  the  baseline,  case  25.  Case  29  shows  a 
background  that  is  highly  time  correlated  but  is  not  related 


TABLE  XIII 


Temporal  and  Spatial  Background  Noise  in  the  Truth 
Model  {S/N  =  2):  Case  26,  27,  28  and  29  Results 


Case  26  Case  27 

.95  Time  Cor.  .5  Time  Cor. 
Spatial  Cor.  Spatial  Cor. 


Average  Standard  Deviation 

X 

y 

x  y 

Atmospheric  Jitter 

Actual 

.22 

.20 

.20  .20 

Filter-Indicated 

.19 

.19 

.19  .19 

Target  Position 

Actual 

19.0  12 

.0 

12.0  9.0 

Filter- Indicated 

.34 

.34 

.34  .34 

Average  Absolute  Error 

Atmospheric  Jitter 

.06 

.05 

.05  .05 

Target  Position 

1.0  1 

.0 

1.0  1.0 

Loss  of  Track 

12  times 

12  times 

Case 

28 

Case  29 

0  Time 

Cor. 

.95  Time  Cor. 

Spatial 

Cor. 

No  Spatial  Cor 

Average  Standard  Deviation 

X 

y 

x  y 

Atmospheric  Jitter 

Actual 

.20 

.20 

.19  .19 

Filter-Indicated 

.19 

.19 

.19  .19 

Target  Position 

Actual 

2.0  1 

.6 

4.0  1.6 

Filter-Indicated 

.34 

.34 

.34  .34 

Average  Absolute  Error 

Atmospheric  Jitter 

.05 

.05 

.05  .05 

Target  Position 

1.0  1 

.0 

1.0  .50 

Loss  of  Track  1  time  8  times 


f 
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spatially  (an  unrealistic  background) .  Filter  performance 
is  poor  but  better  than  that  displayed  by  cases  26  and  27 
and  worse  than  the  performance  of  case  28.  In  summary,  at 
low  S/N  temporally  and  spatially  correlated  background  noise 
can  significantly  hinder  filter  performance.  These  cases 
suggest  that  the  temporal  noise  correlations  may  be  more 
critical  than  the  spatial  noise  correlations.  This  may  not 
be  true  in  reality,  because  the  spatial  correlations  were 
fairly  realistic  (see  Appendix  B)  whereas  the  temporal  cor¬ 
relations  were  worst  case  and  possibly  far  from  reality  (see 
Appendix  C) .  Target  motion  will  also  tend  to  decrease  the 
importance  of  temporal  correlations. 

The  last  two  cases  tested  were  run  with  S/N  equal  to  2 
with  the  target  moving.  The  motion  scenario  previously 
described  in  this  chapter  and  Chapter  II  was  diagonal  up 
and  to  the  right  with  an  initial  velocity  of  28.2  pixels 
per  second.  Case  30  was  run  without  temporal  and  spatial 
background  correlations;  case  31  was  run  with  them.  The 
results  are  shown  in  Table  XIV.  For  both  cases  the  combina¬ 
tion  of  target  motion  with  low  S/N  resulted  in  almost  equal, 
unacceptable  filter  performance.  The  cases  provided  little 
new  insight  to  the  filter.  With  motions  greater  than  eight 
pixels  per  sample  period  which  is  very  easily  reached  for  most 
air-to-air  missile  scenarios,  temporal  modelling  of  the  back¬ 
ground  changes  completely  in  a  single  sample  period. 

Major  Problem  Areas 

The  results  of  the  computer  simulations  presented  in 
tnis  chapter  have  provided  considerable  insight  to  the  weak¬ 
ness  of  the  original  filter.  From  the  parameter  mismatch 
'ins  the  following  list  of  major  problem  areas  in  order  of 
•port ance  resulted: 

Low  signal  to  noise  ratio 

’•r get  size  mismatch 
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TABLE  XIV 


Temporal  and  Spatial  Background  Noise  with  Target 
Motion  (x  =y  =20  pixels/sec)  in  the  Truth  Model  S/N  =  2 
Case  30  and  31  Results 


Case  30  Case  31 

Baseline 

0  Time  Cor.  .95  Time  Cor. 

No  Spatial  Cor.  Spatial  Cor. 

x  y  x  y 

Average  Standard  Deviation 
Atmospheric  Jitter 

Actual  .19  .19  .19  .19 

Filter-Indicated  .19  .19  .19  .19 

Target  Position 

Actual  7.0  10.0  10.0  10.0 

Filter- Indicated  .34  .34  .34  .34 

Average  Absolute  Error 

Atmospheric  Jitter  .08  .05  .06  .08 

Target  Position  7.0  10.0  6.0  12.0 

Loss  of  Track  17  times  18  times 

3.  Target  intensity  mismatch  (very  much  related  to  the 
low  S/N  problems) 

4.  The  ratio  of  rms  target  dynamics  to  rras  atmospheric 
jitter  mismatch. 

From  the  new  truth  model  cases  another  list  of  significant 
problems  are  noted  as  follows  (in  decreasing  order  of  concern) 

1.  Target  motion 

2.  Low  S/N  ratio 

3.  Target  shape  and  orientation  mismatch 

4.  Temporal  and  spatial  background  noise  at  low  S/N 
(but  at  moderate  expected  S/N  levels  these  effects 
were  not  as  significant  as  others) . 


The  various  solutions  to  these  problems  range  from  sim¬ 
ple  to  very  complex.  The  low  signal  to  noise  problem  has 
been  noted  in  many  systems  dealing  with  electromagnetic 
radiation.  It  was  felt  that  the  scope  of  this  research 
should  be  limited  to  exclude  this  major  problem  area.  The 
problem  areas  in  this  research  which  appeared  significant 
and  could  be  addressed  within  the  allotted  time  period  were 
target  motion,  target  size  and  shape,  and  target  intensity. 
The  solution  of  the  target  intensity  problems  will  help 
solve  the  S/N  ratio  problems  noted  in  this  chapter.  Possi¬ 
ble  filter  changes  to  maintain  track  on  a  moving  target 
cover  a  large  range  of  models  and  techniques,  both  adaptive 
and  ad  hoc,  and  will  encompass  much  of  the  rest  of  this 
research.  For  the  problems  caused  by  target  size,  shape, 
and  intensity,  the  apparent  solution  is  real  time  adaptive 
estimation  of  the  target  parameters.  Adaptive  estimation 
schemes  for  estimating  these  parameters  range  from  such  com¬ 
plicated  formulations  as  the  solution  of  the  full  maximum 
likelihood  equations  to  simple  approximations  and  ad  hoc 
schemes  (Ref  16:10).  These  techniques  will  be  studied  in 
this  research,  attempting  to  achieve  a  final  design  with 
desirable  performance  and  yet  online  applicability. 
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IV.  Filter  Requirements 


Introduction 

Section  III  was  designed  to  indicate  the  changes  needed 
to  Capt.  Mercier's  filter  to  enable  it  to  track  near  range 
moving  targets.  Several  areas  needing  attention  are  as 
follows:  (1)  The  filter  must  have  the  ability  to  predict 

future  (next  sample  time)  position,  as  the  missile  can  travel 
many  times  the  width  of  the  FOV  in  one  sample  period.  At 
10  km  a  typical  missile  can  travel  approximately  80  pixels 
per  sample  period.  With  an  8  x  8  FOV  and  no  prediction 
capability,  the  target  image  will  depart  the  FOV,  resulting 
in  loss  of  track.  Therefore,  the  tracker  must  have  a  pre¬ 
dicted  position  for  the  next  sample  time.  (2)  The  filter 
must  be  able  to  estimate  orientation,  size  and  shape  of  the 
image.  This  is  necessary  because  no  a  priori  estimate  of 
flight  path  is  possible  and  also  missiles  differ  in  length 
and  aspect  ratio  (lengthrdiameter) .  (3)  The  filter  must 

determine  the  extent  to  which  a  particular  flight  path  dif¬ 
fers  from  the  dynamics  model  in  the  filter,  in  other  words, 
establish  (i.e.  adaptively  estimate)  the  strength  of  the 
driving  noise  for  its  state  stochastic  process  model.  Flight 
path  identification  can  be  viewed  as  essentially  maneuver 
detection  followed  by  adjustment  of  the  noise  strength  to 
fit  the  trajectory,  and  alternative  means  of  responding  to 
detected  maneuvers  should  be  considered.  (4)  The  filter 
must  have  a  good  estimate  of  the  maximum  intensity  of  the 
image.  As  the  residuals  are  the  key  to  any  parameter  iden¬ 
tification  technique,  a  large  difference  between  the  filter's 
assumed  I  (Ir)  and  the  true  I__„  will  yield  residuals 
with  inaccurate  information.  Each  of  the  above  mentioned 
requirements  and  the  basic  Kalman  filter  equations  will  be 
discussed  individually  in  the  following  paragraphs. 

For  the  filter  to  predict  future  position,  a  current 
estimate  of  velocity  is  required  as  a  minimum.  Using  the 


predicted  position,  the  controller  was  assumed  to  be  able 
to  point  without  error.  Actually,  a  filter  without  predic¬ 
tion  capability  could  track  a  moving  target  if  it  did  not 
leave  the  FOV  in  a  sample  period.  However,  the  pointing 
command  would  still  be  one  sample  behind  at  all  times.  At 
best,  without  predicting’ future  position,  the  tracker  would 
lag  the  missile  considerably  and  in  many  cases  lose  it 
entirely. 

The  six  state  filter  will  be  described  first  and  in 
some  detail.  The  eight  state  filter  will  be  presented 
briefly  showing  the  differences  between  the  two  filters. 


Six  State  Filter 

The  dynamic  model  for  the  six  state  filter  is 


*D<t) 

it 

5 

(94) 

Yo(t) 

n 

ft 

(95) 

where  the 

subscript  D  stands 

for  dynamics,  and 

EtWpCt)] 

=  0 

(96) 

E [Wp  (t)Wg(s)J 

=  a*  6 (t-s) 

(97) 

Equations  (94)  and  (95)  above  are  not  strictly  correct,  as 

the  tracker  orientation  is  fixed  with  the  LOS  which  is  a 

rotating  reference  frame.  However,  a  sample  period  of  1/30 
second  was  considered  short  compared  to  any  time  constant 

describing  missile  acceleration.  In  other  words,  inertial 
equations  were  considered  adequate  for  one  sample  period. 

To  correct  the  filter  model,  equations  similar  to  (63)  and 
(65)  in  the  truth  model  section,  or,  at  least  Coriolis 
theorem  terms  would  need  to  be  incorporated.  Explicit  cor¬ 
rection  of  the  rotating  frame  effects  would  only  be  incor¬ 
porated  if  warranted  in  a  performance/computational  loading 
tradeoff. 
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The  atmospheric  disturbance  model  for  each  axis  in  the 
filter  is 

xA(t)  =  xA(t)  +  wA(t)  (98) 

A 

where  the  subscript  A  stands  for  atmospherics,  and 

ta  =  correlation  time 

E[wA(t)]  =  0  (99) 

2c2 

E[wA(t)wA(s)3  =  6  <t-«)  (100) 

A 

For  a  detailed  explanation  of  this  model  see  Ref  5  and 
Appendix  A. 

Combining  the  dynamic  and  atmospheric  models  yields  a 
propagation  equation  upon  which  to  base  a  Kalman  filter  of 

Xp  *  F(t)Xp(t)  +  Gp  ( t )  Wp  ( t )  (101) 

where  the  subscript  F  stands  for  filter,  and 
F (t)  «  F  *  filter  plant  matrix 
Xp(t)  =  filter  state  vector 
x^  =  azimuth  position 
x2  *  elevation  position 
x3  ®  azimuth  velocity  *  x^ 
x^  *»  elevation  velocity  =  Xj 
Xj  *  azimuth  atmospheric  disturbance 
Xg  *  elevation  atmospheric  disturbance 
Gp(t)  *  Gp  =  6  x  4  noise  distribution  matrix 
E[Wp(t)]  »  0 


71 


ElWpftJWpte) J 


Qp  5( t-s) 


0 

0 


ta 

0 


2cJ 


fi(t-s) 


(102) 


2 

aD  =  strength  of  velocity  noise 
2 

oA  *  variance  of  atmospheric  jitter 
model  output 


The  velocity  noise  strength,  a£,  need  not  be  equal  for  the 
azimuth  and  elevation  axes,  conceptually.  For  this  research, 
the  axis  noise  strengths  were  set  equal  for  the  nonadaptive 
phases  of  tracking  (acquisition)  and  allowed  to  be  unequal 
during  adaptive  phases  of  tracking. 

Explicitly,  equation  (101)  is 


x 


0 

0 

0 

0 

0 

0 


1 

0 

0 

0 

0 

0 


0 

1 

0 

0 

0 

0 


0 

0 

0 

0 

-1 


0 

0 

0 

0 

0 

-1 


0 

0 

0 

1 

0 

0 


0 

0 

0 

0 

1 

0 


w. 


w- 


W, 


w, 


(103) 


where  w^  through  w^  are  independent  white  Gaussian  noises. 
The  filter  state  and  covariance  propagation  equations  that 
result  from  this  propagation  equation  are 


~*fci+l) 


♦p (ti+i,ti)x(ti) 


(104) 
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with  covariance  of 


-^i+1*  “  “F^ti+l,ti^-(ti)^F(ti+l,ti)  + 


j  1  1  £p(ti+1,T)GF(T)QF(T)G^(T)^(ti+1,T)dT  (105) 


The  -  and  +  superscripts  refer  t.o  before  and  after  a  meas¬ 
urement  update  is  made.  Continuing,  the  state  transition 
matrix  is 


Note  that  is  a  constant  matrix.  The  integral  term  in 
equation  (105)  is  evaluated  and  becomes  a  representation  of 
the  growth  in  uncertainty  due  to  driving  noise  since  the 
last  measurement  update.  This  integral,  also  a  constant 
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so  the  covariance  propagation  equation  becomes 


=  +  Spa  (108) 

The  above  derivation  can  be  verified  with  any  text  on  Kalman 
filtering  (Ref  8:171). 

Eight  State  Filter 

The  six  state  filter  above  is,  to  be  sure,  not  the  only 
possibility  for  filter  formulation.  An  eight  state  filter 
including  acceleration  in  both  directions  would  seem  to  be 
appropriate  also.  For  eight  states  the  dynamic  model  would 
be ,  for  example 


=  Hd 

(109) 

• 

2b 

“  2d 

(110) 

• 

2d 

=  «d 

(111) 

where  w_  is  a  vector  of  zero  mean  white  Gaussian  noise. 

— D 

Another  model  which  allows  for  shorter  time  correlation  of 
the  acceleration  is 

a  =  ^  a  +  w_.  (112) 

—  T  —  — D 


where  is  a  vector  of  zero  mean  white  Gaussian  noise  and 
r  is  the  correlation  time.  This  model  will  give  more  rapid 
variations  in  a  than  equation  (111)  but  introduces  another 
parameter,  t,  to  estimate. 

A  third  possible  model  for  a  is 


where 


-o>  v  +  Wp 


v  =  velocity 
|v  x  a  | 


(113) 


b) 


V 


(114) 


and  Wp  is  a  vector  of  zero  mean  white  Gaussian  noise.  This 
model  is  more  complex  and  computationally  less  attractive, 
but  it  forms  the  basis  of  a  model  suggested  for  enhancing 
close  range  tracking  (Ref  17). 

Consideration  was  given  to  these  models  and  equation 
(111)  was  chosen.  Equation  (114)  is  better  suited  to  an 
aircraft  rather  than  a  missile.  Combined  with  increased 
computation  load,  it  seemed  less  appropriate  them  the  other 
models.  With’  adaptive  estimation  of  the  driving  noise, 
equation  (111)  would  give  adequate  performance  in  low  and 
high  acceleration  cases  without  the  need  for  equation  (112) . 

Some  advantages  of  the  eight  state  filter  are  better 
characterization  of  the  true  acceleration  than  the  previous 
(six  state)  model  and  better  prediction  capability.  The 
price  paid  for  this  is  longer  computation  time  and  more 
memory.  Other  advantages,  associated  with  noninertial 
acceleration  compensation  and  maneuver  detection,  will  be 
discussed  subsequently. 

The  statistics  of  the  dynamic  driving  noise  are 

ElWgtt)]  =  0  (115) 

Elw^tJw^CsJJ  =  °D  —  5tt"s)  (116) 

where  the  noise  now  enters  at  the  jerk  (time  derivative  of 
acceleration)  level.  Again,  this  model  is  not  entirely  cor¬ 
rect  either,  but  is  a  better  approximation  than  that  of  the 
six  state  filter.  The  development  of  the  eight  state  filter 
parallels  the  development  of  the  six  state  filter  so  only 
the  results  will  be  presented  here.  The  state  vector  is 

x ^  =  azimuth  position 

Xj  *  elevation  position 

x^  «  azimuth  velocity  =  x^ 

x^  =  elevation  velocity  =  x2 


The  advantages  and  disadvantages  of  both  six  and  eight 
state  formulations  were  considered  and,  unless  performance 
required  the  eight  state  filter  to  be  used,  the  six  state 
filter  would  be  pursued.  Thus,  the  six  state  filter  above 
is  the  basic  formulation.  The  adaptive  estimation  of  the 
QFd  matrix  and  other  modifications  of  this  basic  structure 
will  be  discussed  later  in  this  section. 

Size  and  Shape  Estimation 

The  second  major  requirement  is  to  determine  the  size, 
shape,  and  orientation  of  the  image  ellipse.  Figure  19 
shows  possible  variables  to  describe  the  image.  The  sub¬ 
scripts  v  and  pv  refer  to  the  image  plane  velocity  direction 
and  perpendicular  to  that  velocity  direction,  respectively. 


Figure  19.  Image  Geometry  Characteristics 


Define 


AR  =  aspect  ratio 


The  measurement  update  for  this  Kalman  filter  is  so 
intimately  related  to  the  estimation  of  cv,  apv  and/or  AR 
that  they  will  be  treated  together. 

The  assumed  measurement  relation  for  an  Extended  Kalman 
Filter  (EKF)  is 


z/t^  =  htxtt^,^)  +  v(ti)  (120) 

where 

Etv/t;^] 

Elv(ti)vT(tj)] 

and  h  is  a  nonlinear  function  relating  the  states  to  the 
measurements.  Specifically,  for  the  jth  of  64  pixels 


=  0  (121) 

-  R(ti)«ij  (122) 


hj ,t±) 


PD  A 


lF(ti)expJ-ilAxviyv) 


r°l 


P3 


[AxvAyv] 


(123) 


where 

=  area  of  pixel  j,  j=l,  ...  64 

Ip  =  (filter  estimate  of)  maximum  intensity 

Axv»  AYV  have  similar  meaning  as  in  eqns.  (30)  and 
(31) ,  and  are  measured  with  respect  to  the  principal 
axes  of  the  ellipse  in  Figure  19. 

ay  =  (filter  estimate  of)  semimajor  axis  of  ellipse 
in  Figure  19 

opv  =  (filter  estimate  of)  semiminor  axis  of  ellipse 
in  Figure  19 

Equation  (123)  is  evaluated  like  equation  (7)  in  the  truth 
model  section  except  that  instead  of  averaging  evaluations 
of  the  integrand  over  16  interior  points  of  each  pixel, 
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only  the  midpoint  is  used.  This  is  done  to  reduce  the  com¬ 
putational  loading  in  the  filter.  Also,  since  the  angle  0 
in  Figure  6  is  not  available  to  the  filter,  it  is  never 
computed.  Instead  cos 6  and  sin0  are  estimated  from  x^ 
(azimuth  velocity)  and  x^  (elevation  velocity)  under  the 
assumption  that  the  major  axis  of  the  ellipse  is  colinear 
with  the  image  plane  velocity  (v^LOS  of  Figure  7) . 


The  measurement  update  relations  for  an  Extended  Kalman 
Filter  using  the  inverse  covariance  form  are 


P“1(tt)  =  P_1  (tT)  +  HT(ti)R-1(ti)H(tjL) 


-1,.  +  ,  ,-l 


p  (t : )  -  [p  A(tTn 


—  i 


K(t±)  =  P(t^)HT(ti)R_1(ti) 


(124) 

(125) 

(126) 


*  , .  + 


x(t±)  =  x(ti)  +  K(ti){£(ti)  -  h(x(t.)  ,  ti) }  (127) 


where  £(t. )  is  a  sample  of  z  at  time  t. . 


The  R  matrix  is  the  covariance  of  the  measurement  model. 
For  pixels  with  assumed  independent,  equal  variance  noise, 
R  will  be  constant  in  time  and  of  the  form 


R  =  Rpl_  and  R 
The  H  matrix  is  defined  as 


-1 


*F  1 


H(ti)  = 


3h  (x,t) 


ax 


X  = 


xCt±) 


(128) 


(129) 


For  details  of  how  this  derivative  is  computed,  see  Ref  5 
and  the  program  listing  (subroutine  MEASF)  found  in  Appen¬ 
dix  F. 

The  inverse  covariance  form  of  the  measurement  update 
equations  was  used  for  computational  efficiency.  This  for¬ 
mulation  requires  inversions  of  two  matrices  both  of  the 
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state  vector  dimension.  Had  the  standard  Kalman  filter 
equations  been  used,  only  one  inversion  would  have  been 
required.  But,  that  one  inversion  would  be  for  a  square 
matrix  of  the  order  of  the  measurement  vector  dimension,  a 
64  x  64  matrix. 


Returning  to  the  problem  of  identifying  av,  o^v  and/or 
AR,  there  are  several  methods  available.  The  uncertain 
parameters  can  be  included  in  the  state  vector  and,  if 
modelled  as  Brownian  (as  in  this  research) ,  requires  statis¬ 
tical  knowledge  of  the  uncertain  parameters,  which  is  often 
difficult  to  provide.  Also,  the  appropriate  driving  noise 
strength  for  the  dynamics  model  may  be  difficult  to  inter¬ 
pret  physically . 

Complete  maximum  likelihood  equations  generally  yield 
good  results  but  could  never  be  implemented  on  line  (Ref  16: 
10-126).  Multiple  model  adaptive  estimation  (Ref  16:10-103) 
is  a  possibility,  as  are  some  ad  hoc  procedures  (Ref  16:10- 
121).  Multiple  models,  however,  would  require  many  models 
to  span  the  variation  of  missile  size  with  range  and  pro¬ 
jected  view  with  respect  to  LOS.  The  method  finally  chosen 
was  a  gradient  minimization  technique.  It  approximates 
full  scale  Maximum  Likelihood  estimation  of  states  and 
parameters,  but  with  significant  computational  advantage, 
critical  to  online  implementation  (Ref  10:10-1). 


Equation  (123),  the  filter  measurement  relation,  can 
be  rewritten  with  the  parameters  of  interest  in  evaluated 
form  as 


!»<2L(ti>'<VAR,e,ti)  "  IrexP 


'( (x-x  ) cos 8+(y-y  )  sine) 


(130) 
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Note  that  only  two  of  cv,  cpv  or  AR  need  be  estimated  as 
they  are  related  by  equation  (119).  For  simulation  pur¬ 
poses,  AR  and  ov  were  estimated  so  as  to  be  directly  com¬ 
pared  to  similar  parameters  already  in  the  truth  model. 

For  actual  implementation,  a  tradeoff  analysis  between  per¬ 
formance  and  efficiency  would  have  to  be  made  to  determine 
which  of  the  three  possibilities  (ov  and  AR,  ov  and  cpv, 

a  „  and  AR)  is  best, 
pv 

The  approach  used  to  estimate  ov  and  AR  was  a  gradient 
recursion,  minimizing  the  quadratic  cost  function  given 
below  in  equation  (131) . 

L(x(t^)  ,ov,AR, 9 , t^) 

*  ir(t.)  -  h(x(t?)  ,t.)]T[r(t.)  -  h(x(tT)  ,t. ) ]  (131) 

JL  X  i  i  X  X 

This  equation  can  be  viewed  as  a  simplified  form  of  maximiz¬ 
ing 

L’dK^J.A)  =  In  fZ(ti)|  Zlt^)  ,  A(-(ti}  (131a) 

where  f  is  the  density  function  of  the  current  measurement 
based  on  all  previous  measurements  and  the  parameter  vector 
(A) .  Explicitly  equation  (131a)  is 

L'  trCt^  ,A)  -  ln[2ir"in/2|PL(ti)  |“1/2] 

-  |ir(ti)-h(x(tT),ti)]TP£1(ti)[r(ti)-h(x(t“),ti)]  (131b) 
where  P^t^  *  H(ti)P(tJ)HT(ti)  +  R(t^) 

The  conditional  density  function  of  equation  (131a)  will  be 
maximized  when  equation  (131b)  is  maximized.  The  first  term 
in  (131b)  does  not  affect  the  maximization  and  in  the  second 
term  PT1  is  deleted  primarily  because  of  low  confidence  in 
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it.  Also,  having  64  measurements  may  tend  to  smooth  the 
function.  Removing  p”^  was  assumed  not  to  change  the  posi¬ 
tion  of  the  minimum  substantially.  Then  the  estimator 
equations  are 


ARp(t.) 

*  ARF(ti_1) 

♦«& 

(132) 

AR  =  ARj, 

SvF(ti) 

^VF(ti-l) 

V 

(133) 

where 

c  =  a  constant  determined  empirically  which  is  the 
same  in  (132)  and  (133) 

Another  advantage  of  this  method  is  that  part  of  the  compu¬ 
tation  needed  for  the  partial  derivatives  in  equations  (132) 
and  (133)  is  already  done  in  evaluation  of  equation  (129) . 

All  three  partial  derivatives  (-^  ,  ,  ~~)  are  computed 

in  the  MEASF  subroutine.  See  the  program  listing  for 
details  of  equations  (132)  and  (133) . 

Although  AR  and  av  were  estimated  every  sample  period 
here,  this  may  be  avoided.  AR  and  ctv  do  not  change  rapidly 
and  could  be  estimated  at  a  slower  rate,  for  example  every 
other  sample  or  every  five  samples.  One  caution  should  be 
noted,  however.  Other  portions  of  the  filter  algorithm 
estimate  various  parameters  and  these  estimations  are  also 
dependent  on  the  residuals  given  by 

r(t±)  *  nt^  -hCx(t7),t.,)  (134) 

The  residuals  are  in  turn  dependent  on  h(x(t^) ,t^)  and 
therefore  AR  and  cv.  So,  equations  (132)  and  (133)  must 
acquire  correct  values  rapidly  and  maintain  them.  If  AR 
and  oy  estimates  differ  considerably  from  the  true  values, 
the  residuals  may  be  in  error,  making  the  estimation  of  other 
uncertain  parameters  difficult. 


/ 


Dynamic  Driving  Noise  Estimation 

The  performance  of  any  Kalman  filter  depends  on  ade¬ 
quate  knowledge  of  certain  characteristics.  Some  of  these 
that  must  be  known  are  the  statistics  describing  the  dynamic 
driving  noise.  Several  methods  to  estimate  the  driving 
noise  adaptively  are  maximum  likelihood,  multiple  model 
adaptive  (Bayesian)  techniques,  correlation,  and  covariance 
matching  techniques  (Ref  16:Chapter  10).  The  ability  to 
change  in  real  time  was  deemed  necessary  because  an 
air-to-air  missile  can  exhibit  a  wide  range  of  dynamic 
characteristics.  Therefore,  maximum  tracking  accuracy  could 
not  be  achieved  with  a  single  choice  of  QF(3.  Even  with 
adaptive  estimation  of  QF(j»  a  compromise  between  statistic¬ 
ally  accurate  estimates  and  rapid  response  to  a  change  in 
Qpd  is  required.  The  two  concepts  tend  to  work  against  each 
other  and  it  is  therefore  difficult  to  design  an  estimator 
with  both  attributes. 

Associated  with  an  adaptive  QFd  estimation,  there  seems 
to  be  two  logical  divisions  of  state  estimation,  an  acquisi¬ 
tion  mode  and  a  tracking  mode.  To  start  the  filter  algo¬ 
rithm,  an  estimate  of  the  state  vector  and  covariance  must 
be  provided.  Because  of  the  distances  an  air-to-air  missile 
can  travel  in  one  sample  time  (80  pixels  or  so)  good  esti¬ 
mates  of  position  and  velocity  are  required  by  the  filter. 

To  allow  for  errors  in  values  handed  off  to  the  filter,  from 
a  radar  for  example,  an  acquisition  mode  was  used. 

The  acquisition  mode  consisted  of  an  initial  filter 
state  covariance  matrix  P^  and  driving  noise  Qpd/  with  sig¬ 
nificantly  increased  values  over  what  normally  would  be  used 
for  the  tracking  mode.  The  matrix  would  be  a  function  of 
the  error  statistics  of  the  hand  off  device  and  would  be 
a  function  of  known  deficiencies  in  the  dynamic  model  of  the 
filter.  P  was  then  allowed  to  propagate  and  update  normally, 
using  values  appropriate  for  tracking,  thereby  reducing 
the  covariance  matrix  to  its  proper  values. 
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The  Qpd  matrix  in  equation  (107)  is  the  result  of  a 
stationary  process  and  is  normally  a  constant  matrix.  To 
obtain  adequate  performance  for  a  maneuvering  target 
needs  to  be  "larger"  than  that  necessary  for  a  benign  tar¬ 
get.  This  larger  value  for  QFd  will,  however,  degrade  the 
filter  performance  with  the  benign  target.  For  the  best 
performance  possible  in  both  situations,  Qpd  needs  to  be 
estimated  adaptively.  Adaptive  estimation  was  not  done 
during  the  acquisition  phase  (first  .5  sec)  unless  the  noise 

strength  that  would  represent  the  maneuver  exceeded  the  fil- 

2 

ter  noise  strength.  Instead  oD  in  equation  (107)  started  at 

a  value  larger  than  normal  tracking  would  require  (typically 
2  3 

oD  =  600  pixels/sec  for  the  six  state  filter  or  600  pixels/ 
sec5  for  the  eight  state  filter) ,  then  was  reduced  linearly 
to  a  value  characteristic  of  a  benign  target,  at  .5  sec 
(typically  oD  *  150  pixels/sec  or  pixels/sec  depending  on 
the  filter  used) .  This  reduces  the  number  of  parameters 
needing  estimation  immediately  and  precludes  the  possibility 
of  induced  stability  problems.  Also  note  that  reducing  o* 
has  no  effect  on  the  atmospheric  portion  of  Qpd  during 
acquisition. 

After  the  acquisition  phase,  adaptive  estimation  of 
gpd  starts.  Again,  several  methods  can  be  used  (Ref  16: 
Chapter  10) .  A  full  scale  estimator  based  on  maximum  like¬ 
lihood  suffers  from  computational  loading  problems.  This 
makes  it  virtually  impossible  to  implement  on  line  without 
simplifying  assumptions .  Multiple  model  adaptive  estima¬ 
tion  techniques  are  possible  but  would  require  more  storage 
and  computing  capability.  The  estimator  finally  selected 
was  an  ad  hoc  procedure  based  on  equation  10-93  of  Ref  10. 

It  can  also  be  interpreted  as  an  approximation  to  a  maximum 
likelihood  estimate. 

Heuristically,  the  estimator  can  be  developed  as  fol¬ 
lows. 
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P(tL) 


2rd(ti-l) 


(135) 


P(t+)  =  P(t^)  -  K(ti)H(ti)P(t?)  (136) 

solving  equation  (136)  for  P(t^)  yields 

P(t^)  »  P(t*)  +  K(ti)H(ti)P(t?)  (137) 

where  the  term  Kft^H  (t^P  (tT)  is  preserved  intact  because 
it  will  be  estimated  on  the  basis  of  observed  residuals. 
Setting  equation  (135)  equal  to  equation  (137)  and  solving 
for  yields 


— Fd 


K(ti)H(ti)P(t~)+P(t^)-*F(ti,ti_1)P(t^_1)ip(tifti_1) 


(138) 


Now  all  that  is  needed  is  an  estimate  of  the  first  term  in 
equation  (138) .  The  state  update  equation  (127)  can  be 
written  as 


x(tt)  -  x (t7)  *  K(ti)r(ti)  =  Ax(t^)  (139) 

An  estimate  of  E[Ax (ti) AxT (t±) )  can  be  had  by  forming 

E ( Ax ( t . ) AxT (t . ) ]  ^  ^  I  [Ax (t . ) AxT (t . ) ]  (140) 

“  1  ~  1  N  j=i-N+l  3  3 


Assuming  the  filter  is  operating  in  steady  state,  Axft^  is 
a  white  sequence  with  covariance  K (t^)H(t^) P (t^)  that  is 
independent  of  x(t^)  (Ref  8:229). 

K(t.  )H(t.  )P(t7)  «  w  r  [Ax(t .)  AxT(t.)J  (141) 

-  i  -  i  -  i  "  j«i-N+l  "  3  3 
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so  an  estimate  of  Qpd(t^)  becomes 

^Fd(ti)  "  ff  *  [Ax(t.)AxT(t.)]  +  P(t+) 

1  j=i-N+l  3  3  ~  x 

-  iF(ti'ti-l)^(ti-l)iF(ti'ti-l)  (142) 

Rather  than  just  average  Ax(ti)AxT(ti> ,  the  entire  relation 
is  averaged  over  the  most  recent  N  samples. 

§Fd(ti)  "  iff  *  {Ax(t.)AxT(t.)  +  P(tT) 
j=i-N+l  3-3  ~  3 

-  MVtj-i^tj-iJ^tj'Vi)}  (143) 

One  advantage  of  equation  (143)  is  that  all  required 
quantities  are  readily  available  from  the  filter  in  the  form 
needed.  Also,  to  reduce  storage  requirements,  a  fading  mem¬ 
ory  approximation  to  a  finite  memory  was  employed.  That  is 

QFd(t.)  «  a  ^(t.^)  +  tl-a)QFdl(t.)  (144) 

where  Sp^^)  is  the  single  term  in  the  summation  of  equa¬ 
tion  (142)  corresponding  to  j=i  and  a  is  a  parameter  which 
essentially  controls  how  long  old  estimates  of  Qpd  are  main¬ 
tained.  If  a=0  this  corresponds  to  using  current  Qpdl  only, 
while  a=l  corresponds  to  ignoring  current  Qpdl  and  using 
only  old  Qpd«  Low  a  implies  little  confidence  in  previous 
estimates,  while  high  a  gives  slow  response  to  actual  Qpd 
changes . 

Typical  values  of  a  are 

.7  <  a  <  1  (145) 

From  partial  investigation  of  the  sensitivity  of  6pd  to  this 
parameter  and  the  desire  for  statistical  confidence,  a  was 
set  at  .8  for  this  research. 


Estimation 


Imax 

Looking  at  equations  (123)  and  (129) ,  also  seen  in 
Chapter  II,  reveals  that  Ip  plays  an  important  role  in  pro¬ 
viding  accurate  residuals.  Because  of  both  uncertainties 
and  changes  in  I  with  time,  no  a  priori  value  is  adequate 
for  all  probable  cases.  Therefore,  on  line  estimation  of 
Ip  is  required.  For  a  detailed  explanation  of  I[t|av  estima¬ 
tion  see  Chapter  V. 

Several  methods  are  possible:  (1)  find  the  largest 
value  in  the  measurement  array  and  set  Ip  equal  to  it; 

(2)  one  dimensional  curve  fitting  to  several  pixels  adjacent 
to  the  pixel  with  the  highest  observed  intensity  and  then 
setting  Ip  equal  to  the  maximum  of  that  curve;  and  (3)  two 
dimensional  surface  fitting  to  several  pixels  surrounding 
the  highest  pixel  and  then  setting  Ip  equal  to  the  maximum 
of  that  surface.  Method  (1)  was  used  here  because  it  was 
easy  to  implement  and  provided  adequate  performance  for 
initial  investigations.  In  equation  form  the  estimator  is 

I„  (t . )  =  max{z . (t . ) 1 1  a  i  £  64}  t  =  t .  (146) 

*  J  1  j  3 

Improved  estimation  performance  can  be  obtained  via  two 
modifications  to  this  method. 

Because  Imav  does  not  vary  rapidly  (over  a  few  samples) , 
a  smoother  estimate  can  be  obtained  by  averaging  over  time. 
This  equation,  similar  to  equation  (144),  is: 

ip(tj)  =  c  Ipftj^)  +  Cl-c)Ip(t..)  (147) 

where  Ip  is  found  by  equation  (146) . 

Even  with  time  averaging  this  is  still  a  biased  esti¬ 
mate  because  I„(t. )  includes  a  noise  component  and  .Gaussian 

X  X 

shape  effects  which  will  be  explained  further  in  Chapter  V. 
An  error  function  was  determined  empirically  of  the  form 


c 


(148) 


This  error  function  should  depend  upon  IF,  apv,  xpeak'  and 


also  but  they  were  not  included  because,  for  the  tests 

was  not  varied  and  the  dependence 


ypeak 

using  equation  (146),  Iroax 
on  Xpgafc  and  y^ea^  can  be  neglected  as  long  as  xpea]C<4  and 
Ypeajc<4  (intensity  peak  remains  in  FOV) .  Thus,  the  error 
function  was  chosen  to  be  a  function  of  only  ov(t^).  See 
the  program  listing  and  Chapter  V  for  more  details. 


i 
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V.  Adaptive  Filter  Performance  and  Modifications 
Introduction 

The  basic  new  filter  design  described  in  Chapter  IV  is 
evaluated  to  determine  its  performance  characteristics. 

Each  facet  of  the  adaptive  filter,  including  the  parameter 
estimation  techniques,  the  acquisition  scheme,  the  six  state 
filter,  and  the  eight  state  filter,  is  tested  and  its  per¬ 
formance  examined.  First,  the  performance  of  the  techniques 
for  I  „  estimation  are  discussed.  Then  o„  and  AR  estimates 
are  compared  to  the  truth  model  values;  some  differences  are 
noted  with  possible  explanations.  Next,  the  acquisition 
strategy  as  formulated  in  Chapter  IV  is  shown  to  be  satis¬ 
factory  for  the  trajectories  studied  and  allows  the  filter 
to  recover  from  limited  hand-off  errors.  Errors  caused  by 
the  non-inertial  properties  of  the  rotating  LOS  reference 
frame  (non-inertial  accelerations)  are  noted  in  cases  in 
which  the  six  state  filter  is  used  and  are  greatly  removed 
by  the  use  of  the  eight  state  filter.  Designed  originally 
to  allow  tracking  of  a  maneuverable  target,  the  QFd  estima¬ 
tion  process  exhibited  stability  problems  for  some  steady 
state  trajectories  and  continual  loss  of  track  for  scenarios 
of  extreme  dynamic  maneuver  initiation.  Modifications  were 
attempted  to  detect  vehicle  maneuvers  and  enhance  filter 
tracking.  These  modifications  were  tested  and  proved  par¬ 
tially  successful.  (The  input  parameters  and  selected  per¬ 
formance  plots  for  all  cases  discussed  in  this  chapter  can 
be  found  in  Appendices  K  and  L,  respectively.) 

I  Estimation 
max  _ 

The  estimation  of  I  uses  one  or  more  of  three  tech- 

max 

niques.  The  first  technique  as  described  in  Chapter  IV 
chooses  IF  to  be  equal  to  the  largest  value  in  the  measure¬ 
ment  array  for  each  sample  period.  The  second  technique, 
also  described  in  Chapter  IV,  was  the  time  averaging  of  the 
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maximum  measurement  as  illustrated  by  equation  (147) .  The 
third  technique  uses  an  error  function  to  remove  a  bias 
caused  by  a  combination  of  target  size  and  shape  and  an 
upward  shift  of  all  simulated  values  which  was  purposely 
programmed  to  avoid  negative  measurement  values.  Physically, 
negative  FLIR  outputs  are  impossible  assuming  no  instrument 
noise.  So,  after  the  measurements  are  generated  in  the 
truth  model  to  reflect  the  target  intensity  function  and 
after  the  simulated  zero  mean  background  noise  is  added  to 
the  measurements,  all  values  are  biased  upwards  so  that  the 
largest  negative  measurement  is  set  to  .1.  This  process 
actually  models  the  physical  effects  of  the  time  vary  ng 
nature  of  the  mean  of  the  background  noise.  This  las'  tech¬ 
nique,  which  was  briefly  mentioned  in  Chapter  IV,  is  devel¬ 
oped  in  more  detail  in  this  section.  Test  data  will  demon¬ 
strate  that  a  combination  of  the  three  techniques  yields 

the  best  estimate  of  I _ . 

max 

The  maximum  measurement  technique  was  used  in  tracking 

a  target,  with  a  *  3  pixels,  a  .  =  1  pixel  and  I  =25 
3  v  c  pv  max 

moving  diagonally  cross  range,  target  trajectory  2  as 

described  in  Chapter  II.  The  estimate  for  I _ „  in  this 

scenario  is  plotted  in  Figure  20.  The  data  plotted  is  the 

I  estimate  at  every  tenth  of  a  second  for  the  first  three 
max 

seconds  of  a  single  simulation  run.  The  time  average  value 
of  I__v  is  24.47,  compared  to  a  true  value  of  25,  with  a 
standard  deviation  of  1.12.  This  compares  roughly  to  the 
standard  deviation  of  the  background  noise,  aN  =  1.414. 

Two  characteristics  of  the  distribution  shown  in  Fig¬ 
ure  20  are  notable.  First,  the  distribution  of  the  points 
demonstrate  that  the  I  estimates  follow  the  added  noise, 

lUuX 

resulting  in  an  undesirably  large  dispersion  of  values. 
Second,  the  average  I  estimate  is  less  than  the  tr 
value.  The  cause  and  solution  for  this  trend  will  be 
cussed  more  in  the  two  paragraphs  to  follow. 
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Figure  20.  I__„  Estimate:  Maximum  Measurement 

IUcLX 

The  second  technique  uses  a  time  average  of  the  maxi¬ 
mum  measurement  value  to  determine  an  I  estimate.  As 

max 

shown  in  equation  (147)  ,  the  current  maximur:  measurement 

is  appropriately  weighted  and  added  to  a  weighted  value  of 

the  previous  estimate  of  I  .  The  purpose  of  this  tech- 

max 

nique  is  to  produce  a  smoother  estimate  of  I  .  Figure  21 

max 

shows  that  a  better  behaved  trend  is  achieved.  This  graph 
uses  the  same  data  values  as  Figure  20,  but  time  averaging 
is  used  weighting  the  current  measurement  maximum  by  .2  and 
the  previous  I_av  estimate  by  .8.  The  average  1  is  24.53 
with  a  standard  deviation  of  .334.  While  the  average  is 
very  close  to  the  average  arrived  at  without  smoothing,  the 


Figure  21.  I _ _  Estimate:  Averaged  Maximum  Measurement 

max 

standard  deviation  has  decreased  by  a  factor  of  3.  Time 
averaging  appears  to  be  a  very  good  technique  to  use. 

For  both  of  the  techniques  just  discussed,  the  average 
estimated  I  was  about  half  of  an  intensity  unit  less  than 
the  true  value.  This  fact  is  attributable  to  the  summation 
of  two  effects.  As  previously  mentioned,  to  avoid  negative 
simulated  measurement  values  (assumed  not  to  be  possible 
with  the  FLIR  sensor  as  instrument  noise  is  assumed  negli¬ 
gible)  ,  all  measurements  are  biased  positively.  This  is 
explained  by  the  truth  model  measurement  scheme.  After  gen¬ 
erating  appropriate  measurement  values  based  on  the  targets 
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intensity  function  with  I  x  equal  to  25,  background  noise 
of  known  strength  is  added  to  each  measurement.  Where  the 
target  intensity  is  small,  this  process  can  result  in  nega¬ 
tive  values  in  the  measurement  array.  To  solve  this  prob¬ 
lem  all  measurements  are  biased  upward  (positive)  to  make 
the  largest  negative  value  equal  to  .1. 

The  second  effect  is  caused  by  a  combination  of  target 
size  and  shape.  Each  measurement  based  on  the  Gaussian 
intensity  function  as  defined  in  equation  (7)  represents 
the  average  target  intensity  over  the  area  of  the  pixel. 
(The  procedure  for  computing  the  pixel  output  due  to  target 
intensity  in  the  truth  model  is  more  completely  detailed  in 
Appendix  D.)  Ignoring  for  a  moment  the  added  background 
noise,  the  average  intensity  value  of  a  given  pixel  will  be 

less  than  I _  even  if  the  maximum  is  found  within  the 

pixel.  This  is  shown  in  one  dimension  in  Figure  22.  If 
the  area  under  the  Gaussian  curve  is  integrated  from  0  to  1 
li.e.  one  pixel)  and  then  divided  by  the  length  of  the 
interval,  the  resulting  average,  Ift  ,  is  less  than  Imay . 


I(x,y) 


Figure  22.  Average  Value  of  Gaussian  Curve 
Over  Interval  0  to  1  (1  Pixel) . 
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If  the  target  is  small  in  size,  i.e.  cpv  =  1  pixel  and 
=  1  pixel,  Iavg  can  be  significantly  smaller  than  ITTiay 
because  variation  of  the  intensity  function  is  large  over 
one  pixel.  For  a  large  target,  the  curvature  of  the  func¬ 
tion  is  more  gradual  and  consequently  the  average  intensity 
is  closer  to  the  maximum  in  the  pixel.  Conceptually,  Iavg 

above  can  be  expressed  as  a  function  of  I  ,  x _ .  ,  y^  ., 

max  peak  •'peak 

av  and  apv.  This  function  could  then  be  solved  for  Imav, 
but  it  would  be  a  very  complex  process.  Instead,  this  func¬ 
tion  can  be  approximated  by  a  polynomial  fit  in  terms  of  av 
only,  as  described  in  the  following  paragraph. 

The  combination  of  the  upward  bias  to  avoid  negative 

noise  values  and  the  downward  bias  due  to  target  intensity 

function  curvature  resulted  in  an  average  I_  in  Figures  20 

max 

and  21  of  about  24.5  for  a  =1  pixel  and  o„  =  3  pixels. 

pv  v  c 

Varying  the  size  of  the  target  in  the  cross  range  trajectory 
results  in  different  averages  for  I  „  as  shown  in  Table  XV. 
In  Table  XV  the  subscripts  F  and  T  correspond  to  filter  and 
truth  model  values,  respectively.  ovp  is  the  estimated 
standard  deviation  of  the  target  in  the  direction  of  the 
target's  velocity  vector.  cv  estimation  and  the  biases 
noted  between  and  average  avp  in  Table  XV  are  discussed 
more  fully  in  the  following  section.  The  data  used  to  estab¬ 
lish  Table  XV  values  for  ayT  -  3.0  was  not  the  same  data  as 
was  used  in  Figures  20,  21,  23,  and  24. 

TABLE  XV 

I-.  Variation  with  cvF  for  Single  Simulations 
of  Cross  Range  Target  Ccfpv=l  pixel) 

Avg  Error 


°VT 

Avg  ovF 

■^max 

Avg  Ip 

^max" 

1.0 

1.57 

25.0 

22.91 

2.09 

3.0 

3.63 

25.0 

24.57 

.43 

5.0 

4,78 

25.0 

25.00 

.0 
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an  error  func- 


From  the  values  of  the  errors  and  cyp, 

tion  can  be  derived  to  compute  a  correction  for  each  I_  , 

max 

estimate.  A  second  order  polynomial  fit  to  the  three  points 
in  Table  XV  (cvp  =  1.57  pixels,  error  *  2.09;  cvp  =  3.63 
pixels,  error  =  .43;  avF  *  4.78  pixels,  error  =  0)  resulted 
in  the  following  error  function 

e(avp)  =  .136  0yp2  -  1.518  0yp  +  4.18  (149) 

A  second  order  polynomial  was  chosen  because  it  was  simple 
to  use  and  the  coefficients  were  determined  exactly  by  the 
three  points.  A  more  precise  functional  form  could  have 
been  derived  which  would  have  required  more  data  points  and 
would  have  resulted  in  a  more  complex  function.  The  I 

max 

estimates  shown  in  Figure  20  were  recomputed  using  the 

error  function  and  are  displayed  in  Figure  23.  In  Figure  23 

the  maximum  measurement  for  each  sample  period  is  corrected 

with  the  error  function  using  the  filter's  estimate  of  0yp 

for  that  sample  period;  time  averaging  is  not  used.  The 

resulting  statistics  of  the  estimation  process  is  an  average 

Ip  of  24.99  with  a  standard  deviation  of  1.16.  The  bias 

has  been  removed  but  the  standard  deviation  is  still  large  - 

ITnav  estimates  are  still  reflecting  the  background  noise. 

The  error  function  was  optimized  for  I _ ,  =  25,  but  would 

max 

provide  some  correction  if  .  were  other  than  25.  How- 

max 

ever,  a  time  shortage  did  not  allow  a  more  detailed  study. 
Further  analysis  could  probably  determine  an  approximating 
error  function  depending  on  oyp  and  Ip  (the  previous 

estimate  of  Iroax)  and  other  parameters  also  for  an  appropri¬ 
ate  range  of  o„  and 

v  max 

The  best  Imay  estimation  procedure  incorporates  the 
time  averaging  to  smooth  the  trend  in  the  estimates  and  the 
error  function  to  correct  for  the  bias  in  the  average  maxi¬ 
mum  measurement.  Shown  in  Figure  24  are  the  same  data 
points  previously  discussed  with  time  averaging  and  error 
function.  The  average  Ip  is  24.99  with  a  standard  deviation 
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Figure  23.  '  Ijnax  Estimate:  Maximum  Measurement  With 
Error  Function 

of  .299.  These  statistics  and  a  visual  examination  of  Fig 
ure  24  reveal  the  success  of  this  combined  estimation  pro¬ 
cess.  The  bias  is  gone;  the  estimated  I  values  are 

max 

grouped  around  the  Ip  =  25  line.  Also,  the  distribution 
of  points  is  well  behaved,  no  longer  tracking  the  total 
noise  amplitude.  This  combined  procedure  results  in  an 
accurate  Imav  estimate  well  worth  the  extra  computations 
needed . 


te:  Averaged  Measurement  With 
ion 


AR  and  ov  Estimation 

The  gradient  minimization  of  the  approximated  likeli¬ 
hood  function  given  in  equation  (131) ,  provided  accurate 
estimation  of  cv  and  AR.  The  step  size  constant  c  intro¬ 
duced  in  equations  (132)  and  (133)  was  determined  empiric¬ 
ally  and  set  equal  to  .001.  Optimization  of  this  parameter 
was  not  accomplished,  but  .001  yields  good  response  and 
damping  characteristics.  The  test  values  of  c  were  between 
.015  and  .0003  inclusive.  The  larger  value  provided  oscil¬ 
latory  characteristics  (underdamped) j  the  estimate  would 
oscillate  above  and  below  the  average  value  (about  .75 
pixels  p-p) .  The  smaller  value  gave  more  consistent  esti¬ 
mates  but  was  slower  responding  (overdamped) .  The  value 
used,  .001,  gave  the  slightly  underdamped  response  desired. 

Figures  25  and  26  are  graphs  of  AR,  ov,  apv  estimates 
(apv  is  implied  by  eq.  (119))  for  a  missile  on  trajectory  2 
(cross  range).  The  truth  model  values  are:  ART  =  5, 
avT  "  5'  apvT  ~  1*  The  resP°nse  from  the  filter  initial 
values  (AR  =  1,  ayF  *  3)  is  rapid  and  produces  values  ade¬ 
quate  for  tracking  within  .1  seconds.  Steady  state  (when 
initial  transients  are  gone)  is  reached  in  less  than  .2 
seconds  and  the  parameters  remain  virtually  unchanged  for 
the  duration  of  the  run. 

Note  the  steady  state  time  averaged  values  of  avF  and 
ARp  differ  slightly  from  the  truth  model  values.  ARp  is 
biased  downward  .93  pixel  with  a  standard  deviation  of  .16 
and  avp  is  also  biased  downward  .17  with  a  standard  devia¬ 
tion  of  .32.  The  standard  deviations  indicate  that  the 
filter  values  are  most  likely  not  equal  to  the  truth  model 
values.  The  cause  for  this  bias  is  probably  threefold. 

First  is  the  different  way  by  which  the  filter  evaluates  the 
h(x(t^),t^)  function  (equation  (123))  and  the  truth  model 
evaluates  the  measurement  intensity  function  (equation  (7)). 
Averaging  16  equally  spaced  points  inside  the  pixel  (truth 
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Time  (sec) 

Figure  25.  Transient  Response  of  av 
and  AR  Estimation 

model)  yields  a  much  better  approximation  to  the  integral 
than  using  just  the  midpoint  (filter).  Second  is  that 
unbiased  errors  are  not  assured  in  this  form  of  the  maximum 
likelihood  estimation  concept  because  several  simplifica¬ 
tions  and  linearizations  are  made  to  the  complete  maximum 
likelihood  equations  to  allow  real  time  application.  Third 
after  the  truth  model  image  is  formed  by  equation  (32) ,  the 
lowest  value  in  the  array  (most  surely  a  negative  number 
due  to  addition  of  samples  of  zero  mean  noise  to  the  pixels 
at  the  edge  of  the  array)  is  set  equal  to  .1.  This  will 


r( 


Figure  26.  Steady  State  Response  of  ay  and  AR 
Estimation 


increase  the  intensity  of  each  pixel.  Since  no  compensation 
is  added  to  the  filter's  estimation  of  this  intensity,  any 
perturbations  to  the  image  caused  by  the  truth  model  simula¬ 
tion  procedure  will  be  reflected  in  the  ovF  and  AR  estima¬ 
tions  of  the  filter. 

The  original  purpose  of  the  cv  and  AR  estimation  was 
to  enhance  state  estimation  by  closely  fitting  the  intensity 
pattern.  As  long  as  this  is  accomplished,  the  exact  values 
estimated  for  cy F  and  ARy  are  of  secondary  importance.  To 
this  end,  an  in-depth  study  to  determine  the  cause  of  the 
difference  was  not  accomplished.  Figure  27  compares  one 
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Figure  27.  Truth  Model  and  Filter  Image 


typical  frame  of  the  simulated  image  to  the  filter's  esti¬ 
mate  of  that  image.  The  close  fit  between  the  two,  espec¬ 
ially  near  the  intensity  peak,  is  evident.  Truth  model 
parameters  for  this  frame  were  =  5.0,  0pvT  =  1,  ART  =  5, 
Imax  “  and  filter  estimated  parameters  were  avF  =  5.15, 
°pvF  “  1*23,  ARp  *  4.18  and  Ip  =  25.33.  Also,  the  truth 
model  image  is  corrupted  with  noise  as  usual. 

Figure  28  is  a  graph  of  avF  and  0pVF  estimates  with 
truth  model  ART  *  3,  opvT  “  1,  =  3.  Figure  29  is  simi¬ 

lar  to  Figure  28  with  truth  model  AR  *  1,  a^v  =  1,  ay  «  1. 

In  the  latter  case  of  Figure  29 ,  the  avF  and  OpVF  are  erratic 
due  to  a  loss  of  resolution  when  and  opv  are  less  than 
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apvF  *  x'  value  was  3 

ovp  *  +  ,  initial  value  was  3 


0 


1 


Time 


2 


Figure  29.  o„  and  a  „  for  AR  *  1  (Truth  Model) 
v  pv 


two.  The  area  of  the  array  where  the  larger  gradients  exist 
simply  are  not  covered  by  enough  pixels  to  provide  smooth 
estimates.  This  erratic  behavior,  although  uneasing,  did 
not  degrade  filter  performance  noticeably.  Time  averaging 
could  be  used  to  smooth  the  estimates  if  required  by  differ¬ 
ent  circumstances. 


Target  Acquisition  Procedures 


The  acquisition  procedures  as  described  in  Chapter  IV 
worked  well  under  nominal  conditions  and  received  limited 
testing  to  determine  their  flexibility.  The  acquisition 


strategy  basically  focuses  on  two  questions:  to  what  value 
should  be  set,  and  how  should  Qp^  be  used  for  acquisi¬ 
tion? 

For  most  cases,  the  acquisition  procedure  was  to  set 
both  and  QF^  large  and  then  to  use  a  linearly  decreasing 
^Fd*  The  values  of  were  large  to  reflect  the  inaccura¬ 
cies  of  the  handoff  values  from  other  tracking  systems. 

The  variances  for  the  position,  velocity,  and  accelerations 

(when  eight  state  filter  was  used)  on  the  main  diagonal 

2  2  2 

were  set  to  25  (pixels  ) ,  2000  (pixels  /sec  ) ,  and  100 

O  A 

(pixels  /sec^) ,  respectively.  The  variance  of  the  atmos¬ 
pheric  jitter  was  kept  at  its  initial  input  value  of  .2 
(pixels  ) .  Except  for  the  main  diagonal ,  contained  only 
zero  entries.  In  simple  terms,  the  high  values  in  P  cause 
the  filter  initially  to  follow  the  measurements  more  closely 

through  an  increase  in  the  gain  matrix.  Also  for  acquisi- 
2 

tion,  oD,  used  in  computing  Qpd  for  the  eight  state  filter 
per  equation  (107)  and  for  the  six  state  filter  per  equa- 

3 

tion  (118) ,  is  set  initially  large  at  600  (pixels/sec  for 

c 

the  six  state  filter,  pixels/sec  for  the  eight  state  fil¬ 
ter)  .  This  value  was  chosen  based  on  observed  values  in 
six  state  filter  simulations.  This  value  was  arbitrarily 
used  with  the  eight  state  filter  and  retained  because  it 
gave  acceptable  performance.  The  value  of  Op  (the  target 
dynamics  portion  of  Qpd)  was  decreased  linearly  until  Op  = 
150 (pixels/sec3  or  pixels/sec5  depending  on  the  filter  used) 
between  times  0 . 0  and  0 . 5  seconds .  The  high  initial  values 
in  reflect  an  uncertainty  in  the  target  dynamics  upon 
acquisition.  At  hand-off  the  missile  could  be  turning, 
linearly  accelerating,  or  in  constant  velocity  straight 
flight,  and  the  tracker  must  be  prepared  to  handle  any  case. 
Because  Qpd  effects  P  through  equation  (108)  and  conse¬ 
quently  the  gain,  the  initially  large  QFd  again  causes  the 
filter  to  strongly  track  the  measurements  though  less 
directly  than  large  PQ.  This  can  be  viewed  as  extending 


the  acquisition  (maintaining  high  P  and  K)  time  interval 
beyond  that  obtained  when  a  high  is  used  with  a  low  con¬ 
stant  Qp^  for  all  time.  The  linear  slope  used  for  QFd  was 
based  on  intention;  better  schedules  are  probably  devisable 
with  more  extensive  study.  Overall,'  these  two  simple  tech¬ 
niques  performed  well  in  cases  33  through  35  in  which  these 
techniques  were  used  in  conjunction  with  the  six  state  fil¬ 
ter,  as  can  be  viewed  in  their  plots  in  Appendix  L.  Because 
these  cases  were  not  directly  designed  to  test  the  acquisi¬ 
tion  process  (i.e.  hand-off  values  of  velocities  and  posi¬ 
tion  to  the  filter  had  no  errors) ,  little  can  be  concluded 
from  them  other  than  the  fact  that  the  acquisition  schemes 
worked . 


One  test  of  the  filter's  acquisition  flexibility,  case 

32,  was  performed  and  demonstrated  good  filter  stability. 

The  test  scenario  was  the  tracking  of  a  target  moving  cross 

range  using  the  eight  state  filter.  The  range  at  closest 

approach  was  10  kilometers.  The  high  scheme  previously 

discussed  was  used,  but  the  linear  schedule  for  Qpd  was  not 

used  because  it  was  linked  in  the  computer  software  to  the 

— Fd  estimati°n  scheme  which  was  deactivated  for  this  case. 

Rather,  Q_,  was  held  constant  throughout  the  runs  with 
2  5 

an  =  600  (pixels/sec  ) .  The  hand-off  values  of  the  inertial 


-508  meters/sec,  yQ  =  -308  meters/sec. 


velocity  were  xQ 
and  zQ  =  8  meters/sec  compared  to  true  values  of  -500,  -300 
and  0  meters/sec.  These  values  were  chosen  to  be  represent¬ 
ative  of  probable  hand-off  errors  from  instruments  such  as 
radars  and  constitute  a  2.4  percent  error.  The  results  of 
case  32  show  the  flexibility  of  these  procedures,  as  can  be 
seen  in  Figures  30  and  31.  Figure  30  shows  the  x  channel 
target  position  mean  error  and  one  sigma  envelope.  The 
plot  shows  an  underdamped  but  stable  response  to  the 
initial  velocity  errors.  Similarly,  Figure  31  displays 
the  x  channel,  target  velocity,  mean  error  and  1  sigma 
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envelope.  The  filters  response  is  again  stable  though 
under damped.  Similarly,  results  were  achieved  for  y  channel 
target  position  and  velocity  errors  (see  Appendix  L) .  The 
target  was  never  lost  in  the  20  simulations  that  constitute 
case  32.  These  results  give  some  confidence  in  the  high  P 
acquisition  technique  and  indicate  that  the  linearly 
decreasing  scheme  is  not  needed  for  successful  target 
acquisition.  However,  much  more  study  is  needed  to  deter¬ 
mine  the  limitations  of  this  acquisition  technique  or  to 
develop  other  techniques.  Further  study  is  desirable  but 
infeasible  in  the  time  allotted.  The  basic  feasibility  and 
performance  of  these  techniques  have  been  demonstrated. 
Another  case  which  warrants  attention  in  resolving  the 
acquisition  strategy  would  rerun  case  32  (hand-off  errors) 

without  increased  P  to  determine  how  critical  it  is  to 

—o 

successful  acquisition. 

Performance  of  Target  Dynamics  Equations 

The  two  formulations  of  target  dynamics  which  resulted 
in  six  state  and  eight  state  filters  were  tested  to  estab¬ 
lish  their  performance  capabilities.  The  six  state  filter 
was  first  tested  against  the  missile  moving  cross  range 
(see  Chapter  II  for  motion  descriptions) ,  and  the  resulting 
performance  plots  indicated  a  serious  problem  attributable 
to  a  non-inertial  acceleration  effect  caused  by  the  rotating 
reference  frame.  A  circular  flight  path  for  the  target  was 
tested  against  the  six  state  filter  to  confirm  that  the 
problem  was  indeed  caused  by  non-inertial  acceleration 
effects  of  the  rotating  reference  frame,  as  will  be  discussed 
in  more  detail  in  the  subsequent  paragraphs  of  this  section. 
The  eight  state  filter,  as  expected,  solved  this  problem  by 
estimating  the  acceleration,  which  compensates  for  the  rotate 
ing  frame  effects. 

Case  33  shows  the  six  state  filter  tracking  the  target 
moving  diagonally  cross  range  starting  towards  the  lower 
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left  hand  corner  of  each  frame.  QFd  estimation  is  active 
during  this  case  without  ad  hoc  Qpd  bounding  (as  will  be 
explained  more  fully  in  the  next  section) .  Figures  32  and 
33  show  the  x  and  y  channel  target  position  error  perform¬ 
ance  plot  for  case  33.  The  projection  into  the  FLIR  image 
plane  of  the  inertially  constant  velocity  vector  of  the 
missile  changes  with  time  as  the  FLIR  system  and  its  refer¬ 
ence  frame  rotates  to  maintain  the  target  in  the  center  of 
the  field  of  view.  Thus,  a  non-inertial  acceleration  is 
created  by  the  rotating  system  and  frame.  This  unmodeled 
non-inertial  acceleration  manifests  itself  in  a  small  but 
definite  diagonal  trend  (positive  slope)  in  the  mean,  and 
mean  plus  and  minus  one  sigma  curves  of  Figures  32  and  33. 

The  mean  value  curve  is  no  longer  centered  about  the  0  line 
but  about  a  line  with  positive  slope  (about  5  to  10  degree 
slope)  which  intersects  the  0  line  about  2.5  seconds,  the 
point  of  minimum  range.  Figures  34  and  35  show  the  x  and  y 
channel  target  velocity  error  performance  plots  for  case  33. 
The  positive  diagonal  trend  is  again  noted,  and  the 
velocity  estimate  diverges  during  the  last  second  of  the 
run  because  the  estimation  does  not  keep  up  when  the 
non-inertial  acceleration  again  becomes  large.  The  accel¬ 
eration  created  by  the  rotating  reference  frame  is  equal  to 
two  times  range  rate  times  angular  velocity  divided  by 
range,  2p«/p  where  w  is  the  angular  velocity  of  the  LOS. 

This  expression  is  the  angular  equivalent  to  the  linear 
translational  term  found  in  an  equation  in  Reference  15, 
page  31.  Thus,  the  non-inertial  acceleration  becomes  insig¬ 
nificant  at  very  long  range,  and  at  minimum  range  for  the 
straight  cross  range  target  because  r  =  0.  (Note  that  u>  is 
at  its  maximum  at  minimum  range.)  The  maximum  absolute 
acceleration  occurs  somewhere  between  these  two  positions 
and  is  a  function  of  the  three  quantities,  p,  p,  and  u. 

To  confirm  the  hypothesis  that  the  unmodeled  non-inertial 
acceleration  was  causing  the  trends  shown  in  case  33,  cases  34 
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and  35  were  run.  The  target  flight  path  as  viewed  from 
above  followed  a  large  circle  with  a  radius  of  20  kilometers 
and  with  the  tracking  system  located  at  the  circle's  center. 
Such  a  flight  path  would  not  result  in  any  non-inertial 
accelerations  because  the  true  target  velocity  vector  would 
have  a  constant  projection  in  the  plane  of  the  FLIR  field  of 
view  rotating  with  the  target.  The  six  state  filter  tracks 
the  target  extremely  well.  However,  the  appearance  of  nega¬ 
tive  values  on  the  main  diagonal  of  P  due  to  numerics  and 
the  Qpg  estimation  caused  some  problems,  and  several  ad  hoc 
fixes  were  tried  as  will  be  discussed  in  more  detail  in  the 
next  section. 

Thus,  cases  34  and  35  represent  the  tracking  of  the 
circular  flight  path  target  with  different  solutions  to  the 
P  problem.  In  case  34  the  main  diagonal  terms  of  the  newly 
computed  Qpdl  (the  second  term  on  the  right  hand  side  of 
equation  (.144))  were  set  to  .1  when  a  negative  value  was 
detected.  The  other  values  in  the  row  and  column  in  which 
the  negative  main  diagonal  term  resided  were  set  to  zero. 
Similarly,  for  case  35  any  negative  values  on  the  main  diag¬ 
onal  of  P  were  set  to  .1  and  the  corresponding  rows  and 
columns  were  zeroed  out.  Figures  36  and  37  show  the  x  and 
y  coordinate  target  position  error  statistics  for  case  34. 
Gone  is  the  diagonal  trend  in  the  statistics  which  had  been 
noted  in  Figures  32  and  33  for  case  33.  Figures  38  and  39 
display  the  x  and  y  channel  target  velocity  errors  for  case 
34,  and  the  non-inertial  acceleration  effects  are  absent, 
too.  Case  35  results  exhibited  the  same  characteristics, 
with  the  disappearance  of  the  trends  seen  previously  in 
case  33. 

The  eight  state  filter,  which  includes  acceleration 
states  for  the  x  and  y  channels,  was  a  logical  solution  to 
the  problem  of  non-inertial  accelerations  of  the  rotating 
reference  frame.  An  altered  dynamics  model  of  the  six  state 
filter  incorporating  appropriate  correction  terms  is  a  more 
direct  approach  to  deal  with  the  problem,  but  it  would 
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require  range  and  range  rate  data  to  work  effectively  and 
would  result  in  a  nonlinear  propagation  formulation  which 
would  be  more  difficult  to  use  than  the  eight  state  filter. 
In  case  36,  the  eight  state  filter  tracks  a  target  follow¬ 
ing  the  cross  range  trajectory  (the  same  trajectory  as  cases 
32  and  33) .  QF^  was  held  constant  throughout  the  simula¬ 
tions  —  a  acquisition  schedule  was  not  used  and  QF{j 
adaptive  estimation  not  incorporated. 

Figures  40  and  41  show  the  x  channel  target  position 
and  velocity  errors,  respectively.  The  non-inertial  accel¬ 
eration  effects  due  to  the  rotating  reference  frame  are  now 
estimated  in  the  filter  so  that  the  diagonal  trend  of  the 
statistics  seen  in  Figures  27  through  30  of  case  33  is  not 
seen  here.  For  a  more  direct  comparison  to  case  36,  cases 
33,  34,  and  36  should  be  rerun  without  QF^  adaptation. 
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X  CHANNEL  OYNRMICS  ERROR  CS/N=12.5) 
Figure  40.  Case  36:  Target  Position  Mean  Error  and 
1  Sigma  Envelope  . 
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X  CHANNEL  VELOCITY  ERROR  CS/N=12.5) 
Figure  41.  Case  36:  Target  Velocity  Mean  Error  and 
1  Sigma  Envelope 
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Qg,^  Estimation  Performance  Analysis 

The  estimation  of  QFd  as  formulated  in  Chapter  IV  met 
very  limited  success  for  the  trajectories  tested  and  even¬ 
tually  became  more  hindrance  than  help  in  terms  of  reaching 
the  research  goals.  Qpd  estimation  was  intended  to  allow 
filter  adaptation  to  changes  in  target  dynamics.  However, 
it  did  not  prove  to  be  responsive  enough  to  prevent  loss  of 
track  in  high  g  turns.  Other  solutions  to  the  high  g  turns 
were  investigated  as  will  be  discussed  in  the  next  section. 
Filter  stability  problems  also  were  caused  by  the  QFd  esti¬ 
mation  technique  for  certain  target  scenarios.  The  first 
cases  of  interest  are  33  and  37  in  which  a  target  moving 
straight  cross  range  is  tracked  by  the  six  state  filter 
with  and  without  QFd  estimation,  respectively.  A  summary  of 
the  filter  performance  for  both  cases  is  found  in  Table  XVI. 
The  estimation  scheme  worked  successfully  in  case  33, 
resulting  in  very  acceptable  results  slightly  larger  than 
those  achieved  without  QFd  estimation.  Because  Qpd  estima¬ 
tion  is  aimed  at  adaptation  to  changing  missile  dynamics, 
it  was  not  expected  to  improve  on  steady  state  (non-accel- 
eratir. g  missile)  filter  performance  unless  there  was  a  sig¬ 
nificant  mismatch  in  non-adaptive  QFd  parameters.  The  non- 
adaptive  QFd  six  state  filter  was  manually  tuned  to  this 
particular  cross  range  trajectory.  Consequently,  its  per¬ 
formance  was  very  good.  Again,  time  constraints  did  not 
allow  an  extensive  study  of  trajectories  with  increased 
motion  (linear  accelerations  and  a  variety  of  turns)  creating 
significant  mismatches  for  the  non-adaptive  QFd  filter. 

As  mentioned  briefly  in  the  previous  section,  the  study 
of  the  unmodeled  non-inertial  acceleration  terms,  which  was 
noted  as  a  performance  degrading  trend  in  case  33,  led  to 
the  study  of  the  circular  flight  path  trajectory  tracked  by 
the  six  state  filter  in  cases  34  and  35.  Then  the  eight 
state  filter  was  tested  against  the  straight  cross  range 
target  in  case  36.  These  cases,  34,  35  and  36,  were  initially 


TABLE  XVI 


Six  State  Filter  Tracking  Cross  Range  Target 
With  and  Without  QFd  Estimation  - 
Cases  33  and  37 


Case  33  Case  37 

w/Qpd  est*  w/o  Qpd  est. 
x  y  x  y 

Average  Standard  Deviation 
Target  Position 

Actual  .22  .19  .20  .17 

Filter-Indicated  .22  .20  .22  .21 

Target  Velocity 

Actual  3.5  3.0  2.5  7.5 

Filter-Indicated  3.0  3.0  5.0  5.0 


Average  Absolute  Error 

Target  Position  .12  .08  .05  .05 

Target  Velocity  1.5  1.5  1.0  1.0 


plagued  with  the  infrequent  appearance  of  negative  values  on 
the  main  diagonal  of  P  due  to  the  estimation.  What 
seems  to  happen  is  that  when  the  dynamics  model  matches  the 
target  trajectory  (i.e.  six  state  filter  with  circular  tra¬ 
jectory  and  eight  state  filter  with  straight  cross  range 
trajectory) ,  values  in  P  become  small  and  susceptible  to 

T 

going  negative.  This  can  be  caused  by  errors  in  using  AxAx 
as  a  representation  of  KHP  in  equation  (138). 

Several  solutions  to  the  P  problem  were  tried  before 
final  ad  hoc  fixes  were  found.  The  initial  theory  was  that 
the  main  diagonal  terms  of  QFdl,  the  second  term  on  the 
right  hand  side  of  equation  (144) ,  should  be  bounded  above 
at  specified  values  and  below  at  zero.  The  lower  bounding 
was  accomplished  by  testing  the  main  diagonal  terms  of  Qp^ 
for  each  sample  period.  If  a  negative  value  was  found,  that 
value  was  set  to  zero  as  well  as  every  value  in  the  row  and 
column  in  which  it  was  found.  For  example,  if  the  third 
term  along  the  main  diagonal  was  found  to  be  negative,  then 
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all  the  values  in  the  third  column  and  all  the  values  in  the 
third  row  were  set  to  zero.  If  the  computation  of  Qpdl 
resulted  in  a  main  diagonal  value  greater  than  a  prespecified 
value,  say  600,  then  that  value  was  set  to  600  and  the  other 
matrix  values  in  the  row  and  column  of  the  main  diagonal 
term  were  set  to  appropriate  proportionate  values  (maintain¬ 
ing  same  correlation  coefficients  through  the  adjustment) . 

For  example,  the  third  main  diagonal  term  is  computed  to  be 
700.  Then  th&t  term  is  set  to  600  and  the  values  of  the 
third  column  and  row  were  multiplied  by  the  square  root  of 
600/700.  However,  these  two  techniques  did  not  totally 
remedy  the  negative  values  on  the  main  diagonal  of  P.  The 
same  lower  bounding  scheme  used  on  QFdl  was  attempted  with 
P.  The  result  was  a  totally  unacceptable  tendency  for  the 
filter  to  lose  track  of  the  target  (go  unstable) .  Finally, 
two  similar  ad  hoc  procedures  were  tested  and  appeared  to 
work.  One  technique  performs  the  same  lower  bounding  pro¬ 
cess  on  QFdl  as  described  above  except  that  a  specified  pos¬ 
itive  value,  such  as  .1,  is  substituted  for  the  main  diagonal 
bound  instead  of  zero.  The  other  technique  is  the  same  as 
just  described  but  performed  on  P.  Lower  bounding  Qpdl  and 
P  at  .1  was  used  successfully  in  cases  34  and  35,  respec¬ 
tively,  as  depicted  previously  in  Figures  31  to  34.  The 
value  of  .1  was  determined  to  be  a  reasonable  lower  bound 
based  on  the  observed  values  of  P  for  a  variety  of  cases. 
Comparable  cases  using  zero  lower  bounding  on  Qp^  and  P  and 
upper  bounding  on  Qpdl  were  not  run  because  usually  the 
resulting  numerical  problems  would  not  allow  the  computer 
simulation  software  to  run  to  completion. 

As  will  be  shown  in  the  next  section,  the  QFd  estimation 
technique  is  not  responsive  enough  to  handle  the  extreme 
maneuvers  of  an  air-to-air  missile,  and  other  techniques 
appear  to  be  more  promising  for  adapting  to  detected  maneu¬ 
ver  initiations.  This  fact,  combined  with  the  problems  just 
discussed,  cause  serious  doubts  as  to  the  use  of  Qpg  estima¬ 
tion,  However,  it  should  be  pointed  out  that  estimation 

124 


iiiiimi  mm  iw 


deserves  a  much  more  thorough  study,  especially  for  scenar¬ 
ios  between  the  benign  straight  line  trajectories  and  vio¬ 
lent  20g  turns,  the  two  extremes  employed  in  this  study. 
Also,  the  robustness  of  a  constant  QFd  needs  to  be  tested. 

Maneuver  Required  Modifications 

The  estimation  of  Qp(j  adaptively,  as  described  in  Chap¬ 
ter  IV,  was  intended  to  allow  tracking  of  a  maneuvering  tar¬ 
get  while  minimizing  position  errors  on  all  trajectories. 

The  first  runs  of  the  eight  state  filter  on  trajectory  3 
(20g  pullup)  were  disappointing.  Several  problems  and  poten 
tial  solutions  surfaced  and  are  discussed  in  the  following 
paragraphs . 

The  FOV  of  the  8x8  array  (.160  yrads  x  160  prads)  was 
too  small  for  highly  maneuvering  targets.  The  20g  lateral 
acceleration  can  cause  the  image  to  leave  the  FOV  in  three 
samples  (.1  second)  if  the  position  update  is  not  accurate. 

A  related  difficulty  is  resolution  (number  of  pixels  per 
FOV  dimension) .  Using  an  image  size  of  5  pixels  by  1  pixel 
tcvT  =  5,  opvT  =1)/  virtually  all  the  position  information 
for  the  apv  direction  is  contained  within  1  to  2  pixels  of 
the  peak  position.  That  is,  the  high  gradients  occur  within 
that  region.  The  20g  acceleration  can  move  the  image  3  to 
4  pixels  laterally  in  two  sample  periods.  Figure  42  shows 
this  idea.  With  this  much  separation  between  the  image  peak 
intensity  position  and  the  filter's  estimate  of  that  posi¬ 
tion,  the  linearization  used  in  the  update 

3h  (x ,  t ) 

H  =  (149a) 

x  =  x(t?) 

is  no  longer  valid.  With  significant  separation,  as  in  Fig¬ 
ure  42,  two  factors  cause  the  H  function  to  be  in  error. 
First,  as  the  peaks  of  the  intensity  image  and  the  filter  h 

A  _ 

function  separate,  the  evaluation  point  (x(t^))  is  incorrect 
and  the  linearization  of  equation  (149a)  is  in  error. 
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Second,  the  tails  of  a  Gaussian  intensity  function  are  nearly 
constant  and  the  partial  derivatives  of  equation  (149a)  are 
approximately  zero.  These  two  effects  combine  to  suppress 
any  update  when  the  situation  of  Figure  42  occurs.  Even  if 
the  image  is  in  the  FOV,  the  tracker  will  probably  lose  the 
target,  which  makes  detection  and  update  imperative  before 
one  sample  is  past. 

The  Qpd  estimation  is  accomplished  in  a  way  such  that 
^Fd  (t^)  is  used  for  the  next  sample  period  to  time  t^+^. 
hus,  for  the  first  sample  after  maneuver  initiation  at  time 

A 

ti»  Qf<j  is  still  based  on  the  (cross  range)  non -maneuvering 
portion  of  the  trajectory.  Therefore,  the  filter  gains  are 
low  and  the  update  is  small  even  though  the  residuals  are 
high.  This  response  was  insufficient  to  allow  continued 
tracking  and  the  image  simply  accelerated  laterally  out  of 
the  FOV.  Moreover,  on  the  benign  portion  of  the  trajectory, 
the  estimate  had  become  small  (as  it  should  with  little 
dynamic  noise)  and  the  fading  memory  smoother  would  not 
allow  a  rapid  increase  in  entries  of  Q„. . 
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This  problem  required  several  parameters  and  some  ele¬ 
ments  of  the  state  vector  to  be  changed  significantly  in  one 
sample  period  to  maintain  track.  For  example,  the  elevation 
acceleration  value,  Xg  in  the  state  vector,  must  be  changed 
from  0  to  approximately  2400  pixels/sec2  when  the  lateral  (6) 
acceleration  is  applied.  '  QF(j  and/or  the  P  matrix  entries 
need  to  be  increased  to  avoid  losing  track  while  changing 
values  in  the  state  vector.  Table  XVII  shows  some  useful 
variables  that  could  indicate  and  quantify  a  maneuver.  The 
quantity  Ax  Ax  is  used  to  indicate  a  maneuver  where 

AxTAx  *  tr(AxAxT)  (150) 

and 

Ax  =  P(tt)HT(ti)R“1Cti)ir(ti)  -  h(x(t“)  ,^>3  (150a) 

Equation  (150)  was  chosen  as  the  maneuver  indicator  because 
AxAx  was  already  available  from  the  Q^,d  estimation  algorithm. 
In  many  applications  the  residuals,  in  one  form  or  another, 
are  used  directly  to  monitor  filter  performance.  But,  in 
this  case  the  residual  vector  consists  of  64  elements  and  to 
use  them  is  computationally  impractical.  A  comparison  of 

m  m 

Ax  Ax,  the  azimuth  and  elevation  elements  of  H  (£-h)  and  true 
errors  for  time  starting  at  3.9  seconds  is  given  in  Table 
XVII.  Since  Ax  Ax  is  used  for  detection,  Ax  would  seem 
appropriate  to  quantify  a  maneuver.  The  difficulty  with 
using  Ax  is  that  equation  (150a)  contains  P(t?)  which  might 
have  been  artificially  changed  (protection  from  negative 
diagonal  elements)  and  possibly  one  or  more  elements  of  Ax 
could  be  zero.  To  avoid  this  possibility  and  the  computa¬ 
tional  burden  of  forming  the  residual  vector  inner  product, 

HT(r-h) ,  a  six  element  vector,  is  used  to  quantify  the 

*“  T 

maneuver.  With  R  being  a  diagonal  matrix,  H  (£-h)  is  already 

computed  in  the  filter  and  requires  no  new  computations. 

Maneuver  initiation  occurs  at  4.0  seconds  and  shows  up  in 

the  filter  at  4.03  seconds. 


127 


TABLE  XVII 
Maneuver  Indicators 


Prospective 

Indicators 


3.90 

3.93 

3.97 

4.0 

maneuver 

initiation 

4.03 

4.07 

4.10 


.1024  -23. 

.0082 
.0063 
. 0376  !  -  6.6 


i  9263 


28.2 

2.8 

-1.5 


743.6 

30.3 

-110.2 

-575.4 

3105.0 

39.0 

-11.4 


Actual 

Performance 

xpeak 
azimuth 
error  6 

a 

^peak 
elevation 
error  6^ 

-.17 

.17 

-.19 

.05 

-.01 

.01 

-.07 

-.19 

-.07 

1.59 

.18 

5.37 

.37 

12.11 

Note  the  increase  m  Ax  Ax  and  the  elevation  position 
T 

element  of  H  (£-h)  at  time  4.03  seconds.  The  sign  of  the 
T 

H  (£-h)  element  tells  the  direction  of  the  maneuver,  in  this 
case  the  acceleration  is  up  Cpositive)  perpendicular  to  the 
velocity  vector.  The  magnitude  of  3105.0  is  related  to  how 
far  apart  the  image  and  filter  function  h  centers  are.  Note 
that  although  the  upward  acceleration  continues,  the  indica¬ 
tion  that  a  maneuver  has  occurred  is  lost  by  4.07  seconds. 
Therefore,  it  is  imperative  that  the  filter  detect  the 
maneuver  immediately  and  perform  the  appropriate  update  at 
4.03  seconds  with  that  knowledge. 

Upon  maneuver  detection  three  actions  should  be  taken. 
At  least  an  approximation  or  crude  estimate  of  acceleration 
that  can  describe  the  maneuver  must  be  provided  to  the 
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filter.  The  filter  gains  must  be  increased  directly  (rather 
them  allow  them  to  build  slowly  by  estimating  Q?d  with 
resulting  P  increase)  to  force  increased  emphasis  of  meas¬ 
ured  data  over  filter  predicted  location,  damp  out  any 
transients  and  account  for  uncertainties  in  maneuver  quanti¬ 
fication.  Also,  the  FOV  should  be  increased  to  avoid  let¬ 
ting  the  image  escape  the  FOV.  Taking  the  above  steps  should 
allow  the  filter  to  track  the  target  after  initiation  of  a 
maneuver.  Expanding  the  FOV  was  not  done  in  this  research 
for  reasons  mentioned  in  Chapter  II. 


Looking  at  Table  XVII  there  is  a  definite  correlation 
between  &a,  6g  and  the  respective  elements  of  H  (£-h)  .  This 
correlation  can  be  exploited  to  quantify  6  and  60.  If  6 

Gt  p  Cl 

and  6^  are  then  assumed  to  have  been  the  result  of  a  discon¬ 
tinuity  in  the  acceleration,  an  estimate  of  this  change  in 
acceleration  can  be  computed. 

For  |fij  *  |6g|  of  Table  XVII,  the  respective  elements 
of  HT(r-h)  differ  considerably.  HT (£-h)  is  given  by 


ht cr-h) 


Bhtxtt^  ,t A) 


9x 


tr(ti) 


Mxlt^)  ,t/)  ] 


x(t7) 


(151) 


which  is  a  strong  nonlinear  function  of  o  ,  a _  and  0  (see 

v  pv 

eqn.  (130)).  Because  a ^  =  5  and  cpvT  =  1  significant  dif¬ 
ferences  will  exist  in  the  associated  |HT(£-h) |  values. 
Notwithstanding,  there  is  a  definite  relationship  between 
the  position  errors,  |6a|  and  |6g|,  and  their  respective 
HT (£-h)  elements . 

An  analytical  solution  for  5Q  and  6^  is  not  possible 
so  an  empirical  approach  was  attempted  to  obtain  this  quan¬ 
titative  relationship.  Two  coordinate  frames  seem  reason¬ 
able  to  approach  this  problem.  The  filter  uses  the  azimuth/ 
elevation  LOS  system,  but  lateral  acceleration  is  expressed 
roost  easily  in  a  missile  oriented  frame  (av,  apy) .  The 
filter  frame  is  desirable  insofar  as  the  result  is  used  in 
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that  frame.  To  this  end,  the  size,  shape  and  orientation 
parameters  of  the  filter  (cvF,  °pvF»  sine,  cose)  are  trans¬ 
formed  into  a&z  and  a^.  This  is  done  using  the  following 
equations . 


o 


az 


CTel 


(152) 


(153) 


where 


e  = 


(154) 


and  subscripts  a2  and  el  stand  for  azimuth  and  elevation, 
respectively . 


The  functional  relationship  g C - )  desired,  for  example 
6p,  is  of  the  form 


-  9'[BT<r-h>,  «.1> 


(155) 


where  a_„  =  image  dispersion  in  the  azimuth  direction 

ael  ™  dispersion  in  the  elevation  direction 

The  function  g ( • )  is  strictly  a  function  of  the  entire 
T 

H  (T-h)  vector.  However,  if  tracking  correctly,  the  depend¬ 
ence  on  all  elements  except  the  one  associated  with  eleva¬ 
tion  (8)  can  be  neglected. 

Because  of  nonlinearities  and  other  considerations, 
plots  of  log|6g|  vs  log] lHT(£-h)]g|  were  made.  Figure  43 

shows  one  such  plot  for  a  case  with  AR_  =  3.  For  values  of 

i  «  T  ^ 

1 6 1  >  .1  the  high  correlation  with  H  (£-h)  is  preserved. 

For  values  of  |6|  <  .1  the  data  is  noisy  and  the  correlation 


1 


i  ) 
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is  lost.  The  function  sought  from  these  graphs  was  of  the 
following  assumed  form 


log | [HT (£-h) J  p |  =  a  log|6g|  +b  (156) 

where  a  and  b  are  the  adjustable  parameters.  In  both  cases 
(AR^  =  5  and  AI^  =  3) ,  a  was  approximately  constant,  within 
the  resolution  of  the  measurements,  and  set  equal  to  .82. 
Figure  44  shows  the  intercepts  (b)  for  the  two  cases  and. 


Figure  44.  Intercept  for  Equation  (156) 


Equation  (156)  then  becomes 
log | lHT(r-h) ] g|  *  .82  log|6B|  -  .46  oel  +  3.9  (157) 

Solving  for  | 6  ^ | 

(  [log  !HT(r-h)]  +  .46  c  -  3.91) 


1 


This  equation  is  applied  to  the  elevation  direction  and  a 
similar  equation  applied  to  the  azimuth  direction  when  a 
maneuver  is  detected  by  Ax  Ax  Cor  H  (£- h)  entries)  surpass¬ 
ing  a  threshold.  The  6's  are  assumed  to  have  resulted  from 
an  acceleration  not  reflected  by  the  acceleration  state 
estimates.  In  equation  form 

6e  =  i(Aae)At2  (159) 

is  solved  for  Aa^  and  the  sign  of  [H  ( £-h) ] ^  incorporated 
to  give 

2  I  6S1  m 

Aa  =  - S~  •  sgn{  [HT  (T-h)  ]  „}  (160) 

B  At2  ~  ° 

where  1 6  ^ (  is  estimated  using  equation  (158).  Assuming  the 

maneuver  was  detected  at  ti#  the  incremental  changes  in 

acceleration  (Aa  ,  Aa0)  are  added  to  the  acceleration  state 
a  c  ^  + 

estimates  and  the  modified  state  vector  set  equal  to  x(ti_1) . 
This  new  state  vector  is  propagated  from  xCtt^)  to  x(t7) 
and  used  as  usual. 

Incorporating  this  algorithm  into  the  filter  and  run¬ 
ning  some  test  cases  showed  promising  results.  The  filter 
would  track  for  several  samples  beyond  the  point  of  loss  of 
lock  in  the  initial  design  without  the  estimated  accelera¬ 
tion  incorporated.  The  acceleration  obtained  from  equation 
(160)  is  not  sufficiently  accurate  to  allow  continued  track¬ 
ing  without  good  updates  in  the  next  few  samples.  Also,  the 
acceleration  is  not  constant  in  the  azimuth,  elevation,  LOS 
coordinate  frame.  Thus,  the  updates  in  the  filter  must  be 
accurate  and  rapid.  Table  XVIII  shows  a  typical  sample 
result  when  the  acceleration  estimate  is  computed  and 
incorporated  into  the  state  vector  but  the  gains  are  left 
unchanged.  Comparing  Table  XVIII  to  Table  XVII  shows  this 
improved  performance. 
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TABLE  XVIII 


Elevation  Acceleration  Estimate  Only 


Time 

6a 

Remarks 

3.9 

.22 

.17 

3.93 

.20 

.03 

3.97 

.37 

-.02 

4.0 

.31 

-.20 

maneuver  initiated 

4.03 

.31 

.23 

mm 

.55 

.68 

.18 

.08 

4.13 

1.69 

.00 

4.17 

2.99 

-.37 

4.20 

4.58 

-.87 

4.23 

6.86 

-1.15 

To  increase  the  filter  gains ,  the  P  matrix  was  redefined 
by  increasing  the  variances  of  the  acceleration  and  velocity 
components.  For  velocity  uncertainty  the  variance  was 
increased  to  100  and  for  acceleration  uncertainty  the  vari¬ 
ance  was  increased  to  300.  For  implementation  these  vari¬ 
ances  would  be  a  function  of  the  accuracy  of  the  maneuver 
detection  algorithm  and  would  be  similar  to  an  acquisition 
mode.  The  performance  of  the  filter  improved  and  it  was 
able  to  track  the  image  for  about  .5  seconds.  Table  XIX 
shows  the  result  when  the  increase  in  elements  of  P  is  added. 
This  sample  corresponds  directly  to  that  of  Table  XVIII 
(same  samples  of  noise  process  used  in  simulation) .  The 
maneuver  detection,  Aa  computation  and  P  increase  function 
properly.  However,  at  4,2  seconds  two  diagonal  elements  of 
the  P  matrix  become  negative  and  the  affected  rows  and  col¬ 
umns  were  set  to  zero.  This  seemed  to  cause  the  algorithm 


TABLE  XIX 


Elevation  Acceleration  Estimate  and  Increased  P  Matrix 


Time 

l 

6a 

*6 

Remarks 

!  3.90 

.22 

.16 

j  3.93 

i 20 

.03 

3.97 

.38 

.03 

4.00 

.31 

-.21 

Maneuver  initiated 

4.03 

.31 

1.59 

First  x  (4.03  )  (unmodified) 

4.03 

-.46 

.2 

New  x  (4.0 3~ )  (modified)  true 
acceleration  was  2448  pixels/ 
sec2  estimated  value  was 

4.07 

-1.68 

-.34 

4.10 

-1.87 

-.35 

4.13 

-2.38 

-.13 

4.17 

-2.18 

.34 

4.20 

-1.59 

.25 

Two  elements  of  P  (4.20+)  <  0. 
Rows  and  columns  set  to  0.  No 
x^  or  Xg  update. 

4.23 

-1.05 

.38 

4.27 

.15 

-.06 

Element  of  P  (4.27+)  <  0.  Row 
and  column  set  to  0.  No  x^ 
update.  ** 

4.30 

1.29 

-.45 

The  <5a  tripped  maneuver  detector 
and  was  corrected 

4.30 

.78 

-.74 

New  x  (4. 30~) 

4.33 

.32 

-1.56 

4.37 

.57 

-1.36 

4.40 

.17 

-1.15 

4.43 

-.17 

-.63 

4.47 

.09 

-.90 

4.50 

.73 

-1.40 

4.53 

1.49 

-1.69 

4.57 

1.47 

-2.33 

4.60 

1.86 

-2.99 

Filter  function  h  and  image  no 
longer  overlap  significantly 
(lost  track) 

to  become  unstable.  Note  the  increased  error  in  the  azimuth 
direction  (acceleration  was  in  the  elevation  direction) . 

This  occurred  because  the  Aa  was  computed  from  a  | 6  |  and 

i  T  i  ”  ^ 

| [H  (£~h)3a|  with  insufficient  correlation.  Even  so,  the 

filter  was  recovering  from  the  erroneous  Aaa  until  the  P 

matrix  becomes  indefinite  at  4.2  seconds.  This  can  be 

remedied  in  part  by  lower  bounding  the  P  diagonal  terms  by  a 

positive  number  (.1  for  this  research).  Continued  work  in 

this  area  is  warranted,  including  the  incorporation  of  FOV 

expansion  and  contraction. 

Final  Filter  Performance 

The  final  filter  algorithm  included  I  ,  AR  and  c  . 

max  v 

estimation,  fixed  matrix,  and  an  acquisition  scheme  all 
incorporated  into  the  eight  state  formulation.  Recall  in 
Chapter  IV  the  desire  to  use  the  six  state  filter  vs.  the 
eight  state  formulation.  Because  of  the  non-inertial  accel¬ 
erations  the  rotating  LOS  coordinate  system  presented,  the 
eight  state  filter  was  selected  to  reduce  these  effects. 

The  Qpd  adaptive  estimation  proved  to  be  inadequate. 

Time  constraints  prevented  a  thorough  investigation  of  the 
causes  for  its  trouble,  but  it  seemed  to  create  more  prob¬ 
lems  than  it  solved.  It  may  be  useful  for  dynamic  environ¬ 
ments  below  maneuver  detection  threshold,  but  its  perform¬ 
ance  benefit  would  need  to  be  compared  to  the  cost  in  com¬ 
puter  time  and  potential  induced  stability  problems  before 
it  could  be  recommended  for  implementation.  So  the  final 
filter  uses  a  Qp^  high  enough  to  allow  acquisition  and  is 
fixed  at  that  value. 

Performance  of  this  filter  on  trajectory  2  (cross  range) 
is  shown  in  several  figures  starting  with  Figure  45.  This 
figure  shows  the  azimuth  mean  position  error  plus  or  minus 
one  sigma.  Comparing  this  to  Figure  32  shows  that  the 
eight  state  filter  exhibits  little  evidence  of  the  non- 
inertial  effects.  The  maximum  mean  error  is  .175  pixels 


136 


HERN  ERROR  4-1  SIGMA 


E  J  > 
a  v  a 

?  I  W 
Df£ 

t-  U  -J 

KxL 


Ot'O 

02  ’0  00*0 

02*0" 

oro-° 

S13XId 

X  CHANNEL  DYNRM1CS  ERROR  (S/N=  12.5) 
Figure  45.  Case  36:  Target  Position  Mean  Error  and 
1  Sigma  Envelope  -  X  Channel  , 
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(3.5  yrads)  and  the  three  sigma  value  for  this  point  is  .84 
pixels  (16.8  yrads).  Figure  46  shows  the  azimuth  mean 
velocity  error  to  less  than  .17%  of  the  velocity  value. 
Generally,  acceleration  magnitude  will  be  nonzero  at  hand- 
off  and  there  will  be  a  transient  before  the  acceleration 
estimate  is  accurate.  This  transient  is  damped  out  by  the 
acquisition  mode  with  its  higher  P  elements  reflecting  the 
uncertainty  at  hand-off. 

Figure  47  shows  the  filter  estimate  of  the  azimuth 
position  standard  deviation  vs.  the  actual  sample  standard 
deviation  from  the  simulation.  Figure  48  shows  the  corres¬ 
ponding  results  for  azimuth  velocity.  Both  plots  show  a 

slight  overestimation  of  the  variance,  but  nevertheless  good 

2 

performance.  The  dynamic  driving  noise  strength  (o^  =  600) 
was  tuned  to  this  trajectory  (2)  and  the  filter  should  per¬ 
form  well.  However,  limited  time  prevented  an  extensive 
investigation  of  suitable  noise  strengths.  A  more  dynamic 
environment  would  make  a  higher  strength  noise  appropriate. 
If  good  performance  required  a  well-matched  Qpd  then  adap¬ 
tive  estimation  of  Qpd  could  be  used  to  advantage.  Figure 
49  through  Figure  52,  for  the  elevation  direction,  are  simi¬ 
lar  to  the  plots  for  azimuth.  The  same  comments  given  above 
are  applicable  to  the  elevation  direction.  Parameter  values 
for  this  run  are  given  in  Appendix  K  case  36. 

For  the  highly  dynamic  environment  of  trajectory  3 
(20g  pullup) ,  partial  results  are  given  in  the  preceding 
sections  of  this  chapter.  The  Monte  Carlo  runs  were  not 
completely  successful,  so  no  corresponding  plots  were  made. 
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X  CHANNEL  VELOCITY  ERROR  (5/N=  12.5) 
Figure  46.  Case  36 j  Target  Velocity  Mean  Error  and 
1  Sigma  Envelope  -  X  Channel 
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Figure  47.  Case  36:  Target  Position  Standard  Deviation, 
Actual  and  Filter-Indicated  -  X  Channel 
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FILTER  V5.  flCTURL  SIOMfl  PLOT  (S/N  =12.5) 
Figure  52.  Case  36:  Target  Velocity  Standard  Deviation, 
Actual  and  Filter-Indicated  -  Y  Channel 
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VI.  Conclusions  and  Recommendations 


Conclusions 

The  scope  of  this  study  allowed  numerous  conclusions, 
particularly  concerning  the  various  filter  techniques  tested. 
Most  of  the  filter  techniques,  such  as  adaptive  schemes  and 
state  space  formulations,  which  are  addressed  in  these  con¬ 
cluding  remarks,  were  formulated  in  detail  in  Chapter  IV 
and  analyzed  in  Chapter  V. 

The  fairly  simple  techniques  developed  for  the  estima¬ 
tion  of  the  maximum  target  intensity,  I  ,  achieved  very 
good  accuracy.  With  a  true  Imax  *  25,  the  estimation  pro¬ 
cess  averaged  24.99  with  a  standard  deviation  of  .299.  The 
adaptive  process  was  specifically  tuned  to  this  level  of 
target  intensity,  but  serious  degradation  is  not  anticipated 
at  other  target  intensities  using  this  estimation  scheme. 

And,  the  small  price  in  computational  load  is  readily  afford¬ 
able  to  identify  this  critical  parameter  accurately. 

Similarly,  the  simple  gradient  schemes  devised  for  esti¬ 
mating  the  size  of  the  target  in  its  longitudinal  axis,  a , 
and  the  aspect  ratio  of  the  target,  AR,  were  sufficiently 
accurate  and  easily  implemented.  These  estimation  schemes 
performed  so  well  that  further  investigations  into  more 
precise  approximations  to  solutions  of  the  full  Maximum 
Likelihood  Equations  were  not  conducted.  For  a  true  AR  of 
5,  the  estimation  scheme  averaged  4.07  with  a  standard 
deviation  of  .16.  The  large  bias  between  the  true  and  esti¬ 
mated  values  is  well  understood,  and  very  accurate  target 
modelling  and  state  estimation  are  still  achieved  as  is 
discussed  in  detail  in  Chapter  V.  For  a  true  of  5,  the 
estimation  of  cv  averaged  4.83  with  a  standard  deviation  of 
.32.  Again,  the  slight  bias  is  understood  and  causes  no 
problems  in  target  modelling. 

The  eight  state  filter  is  preferred  to  the  six  state 
filter  to  provide  more  accurate  tracking  of  the  missile  as 
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cam  be  seen  by  comparing  Figures  32  through  35  (six  state 
filter)  to  Figures  40  and  41  (eight  state  filter)  in  Chapter 
V.  The  eight  state  filter  is  needed  to  accommodate  the 
non-inertial  accelerations  caused  by  the  rotating  FLIR  ref¬ 
erence  frame  as  seen  in  the  straight  line  flight  test  cases. 
Also,  the  acceleration  estimation  scheme  developed  for  high 
g  turns  is  more  easily  used  in  conjunction  with  the  eight 
state  filter  because  x  and  y  channel  accelerations  are  two 
of  its  states. 

The  technique  for  estimating  the  dynamic  driving 
noises,  Qf<j,  was  not  conclusively  determined  to  be  useful. 
Much  was  learned  concerning  procedures  to  combat  apparent 
stability  problems  that  occur  when  using  the  Qp^  estimation 
process.  Zero  lower  bounding  of  P  and  QF(j/  and  upper  bound¬ 
ing  on  0Fd  were  not  effective  solutions  to  the  instability 
of  the  filter  which  appeared  as  negative  values  on  the  main 
diagonal  of  P.  A  brief  study  of  positive  lower  bounding  of 
P  and  Qp^  demonstrated  good  potential  for  solving  the  sta¬ 
bility  problems  noted.  It  was  discovered  that  Qpd  estimation 
was  not  fast  enough  for  hard  maneuvers  but  it  could  be  used 
perhaps  for  intermediate  (medium  maneuvers)  cases  which  were 
not  investigated  in  Chapter  V. 

The  maneuver  detection  and  acceleration  estimation 

techniques  showed  great  promise  in  leading  to  the  ability 

to  track  the  missile  in  high  g  turns  with  a  very  limited 

field  of  view.  The  maneuver  detection  technique  is  based 

T 

on  large  changes  in  the  scalar  Ax  Ax,  where  Ax  equals  the 

T 

gain,  K,  times  the  residuals,  or  on  large  changes  in  H  [£-h] 

+  T  —  1 

(note  K  *  P  H  R  and  R  *  R  I_  here)  .  These  techniques 

successfully  detected  all  the  maneuvers  tested  which  were 

extreme.  The  acceleration  estimate  which  was  based  on  the 

linearized  measurement  matrix  multiplied  by  the  residuals, 

T 

H  (r-h)  proved  capable  of  identifying  the  acceleration  of 
a  high  g  turn  to  within  five  percent  of  its  true  value  in 
one  sample  period.  Immediately  after  maneuver  detection 
two  actions  were  taken:  P+  was  increased  which  immediately 


increased  the  gain  (K  =  P  HR  A) ;  and  an  acceleration  esti¬ 
mate,  as  described  above  and  in  Chapter  V,  was  incorporated 
into  the  filter.  These  procedures  were  effective  for  about 
half  a  second,  but  the  filter  then  lost  track  of  the  tar¬ 
get.  The  loss  of  track  was  not  directly  attributable  to 
the  initial  actions  taken,  but  was  caused  by  the  filter’s 
inability  to  keep  the  gains  high  enough  throughout  the  turn. 
A  procedure,  which  was  not  evaluated  but  could  be  effective 
in  tracking  missiles  in  g  turns,  would  additionally  expand 
the  field  of  view  of  the  filter  by  summing  up  pixel  values 
to  represent  a  single  measurement,  and  is  described  in  the 
recommendations  section  to  follow. 

Finally,  in  the  high  levels  of  S/N  (10  or  greater) 
studied,  the  spatial  and  temporal  correlations  of  the  back¬ 
ground  noise  were  not  found  to  be  significant  in  their 
effects  on  the  filters  performance,  as  was  shown  in  Chapter 
III.  They  were  found  to  have  a  significant  impact  at  S/N 
equal  to  two  or  lower.  The  combination  of  spatially  corre¬ 
lated  noise  and  temporally  correlated  noise  was  found  to 
have  more  effect  than  each  one  taken  separately  as  can  be 
seen  in  the  cases  presented  in  Chapter  III. 

Recommendations 

Many  areas  needing  more  extensive  study  were  bypassed 
in  an  effort  to  progress  the  general  design  of  an  effective 
filter  as  far  as  possible.  Hopefully,  this  study  has  solved 
some  problems  and  provided  useful  insights  towards  a  final 
adaptive  extended  Kalman  filter.  Future  research  on  this 
problem  can  go  in  many  directions. 

More  study  needs  to  go  towards  testing  the  final  filter 
design  to  determine  more  completely  its  capabilities  and  to 
see  if  it  can  track  a  turning  missile  with  minor  modifica¬ 
tions.  In  particular,  investigations  of  target  scenarios 
between  the  straight  missile  flight  path  and  the  20  g  turn 
are  needed.  Also,  the  maneuver  detection  and  acceleration 
estimation  schemes  will  require  more  study. 
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Some  major  modifications  could  include  (but  are  not 
limited  to)  alternate  reference  frames  for  problem  formula¬ 
tion,  alternate  target  acceleration  models,  enlarging  the 
field  of  view  (dynamic  pixel  compression  schemes) ,  and  dual 
parallel  filters  (a  form  of  multiple  model  filtering.  Ref  8: 
10-103  to  10-106) .  Alternate  reference  frames  could  possi¬ 
bly  lead  to  more  accurate  or  more  efficient  filter  formula¬ 
tions.  For  example,  a  vehicle  fixed  reference  frame  filter 
formulation  Could  provide  better  filter  performance  during 
high  g  missile  maneuvers.  Alternate  models  for  target 
acceleration,  such  as  constant  turn  rate  model  or  first 
order  Markov  process ,  could  enhance  filter  performance  as 
compared  to  the  present  model  which  represents  acceleration 
as  Brownian  motion.  The  idea  of  enlarging  the  field  of  view 
is  motivated  by  the  need  to  keep  a  maneuvering  target  in  the 
field  of  view  considered  by  the  filter.  One  such  scheme 
would  sum  up  2-by-2  or  4-by-4  arrays  of  pixels  into  one 
measurement  to  allow  an  expansion  of  the  field  of  view 
without  changing  the  size  of  the  measurement  array  in  the 
filter.  The  dual  filters  idea  may  be  most  useful.  One 
filter,  such  as  designed  here,  could  be  used  to  track  linear 
missile  trajectories,  and  the  other  filter  with  high  gains, 
possibly  constant,  would  track  the  high  g  missile  turns. 
Switching  between  filters  could  be  provided  by  a  maneuver 
detection  scheme  as  studied  in  this  research. 

Another  study  could  concentrate  on  the  estimation  tech¬ 
nique  for  in  order  to  analyze  and  solve  the  stability 
problems  noted  in  this  research.  The  positive  lower  bound¬ 
ing  of  Qpd  and  P  would  be  good  starting  points  for  such  a 
study,  but  a  greater  variety  of  target  scenarios  should  be 
investigated  to  establish  more  universal  lower  bounds  for 
the  diagonal  terms  of  each  matrix.  Further  research  into 
Qpd  estimation  techniques  should  address  its  role  in  target 
acquisition  strategies,  and  a  tradeoff  study  should  be  con¬ 
ducted  to  determine  whether  performance  benefits  warrant  the 
additional  complexity  of  the  adaptation. 
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An  idea  which  may  merit  investigation  is  the  modelling 
of  some  targets  by  two  or  more  related  elliptical  contour 
intensity  distributions.  Targets  such  as  aircraft,  may 
appear  to  be  composed  of  more  than  one  major  source  of 
infrared  (IR)  radiation  when  viewed  by  an  IR  sensor.  For 
example,  aircraft  nose,  wingtips,  engines  and  engine 
exhaust  are  all  possible  sources  of  high  IR  radiation. 

More  effective  tracking  schemes  could  possibly  be  formulated 
.  with  multiple'  intensity  patterns  which  maintain  a  consistent 
geometrical  relationship  to  each  other. 

A  more  long  range  research  topic  could  focus  on  an 
efficient  on  line  implementation  form  of  proposed  filters, 
such  as  square  root  filter  forms  or  U-D  factorization  forms 
(Ref  15:  Chapter  7).  The  sample  rate  for  any  effective 
tracking  system,  such  as  the  30  Hz  rate  currently  used,  ? 
will  have  to  be  relatively  fast.  This  will  require  extreme 
efficiency  and  accuracy  of  the  on  line  formulation  used,  and 
optimization  of  the  algorithm  in  this  regard  should  be 
undertaken . 

The  very  specialized  problem  of  tracking  at  low  S/N 
will  eventually  need  to  be  confronted.  As  shown  in  the 
present  research,  a  part  of  such  future  study  should  include 
the  spatial  and  temporal  nature  of  the  background  noise. 

Thus,  part  of  the  study  would  entail  more  extensive  real 
data  analysis  to  evaluate  the  temporal  and  spatial  correla¬ 
tions  involved  more  precisely  and  to  formulate  new  models 
as  needed.  This  possibly  could  be  accomplished  by  incorpor¬ 
ating  the  temporal  and/or  spatial  noise  concepts  into  the 
filter,  and,  then,  evaluate  filter  performance  against  a 
realistic,  more  complex  truth  model  such  as  stored  test  data. 

Future  filter  designs  to  incorporate  range  and/or  range 
rate  data  may  provide  the  best  approach  to  the  missile  track¬ 
ing  problem.  Such  measurements  would  allow  direct  computa¬ 
tion  of  compensation  terms  for  non-inertial  acceleration 
effects,  and  would  enhance  estimation  of  range -dependent 
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parameters  to  describe  target  size  effects  and  maneuver 
rate  capabilities.  The  price  of  the  extra  sensors  needed 
for  the  range  and  range  rate  data  could  easily  be  outweighed 
by  the  accuracy  and  simplification  gained. 

The  ultimate  goal  would  be  the  formulation  of  a  sub- 
optimal  closed  loop  controller  design.  Incorporating  both 
laser  system  and  tracking  system  dynamics  into  a  near  opti¬ 
mal  feedback  system  (as  via  forced  separation  design,  sepa¬ 
rating  the  generation  of  deterministic  controllers  and 
estimators,  and  then  combining  them)  would  be  a  very  complex 
problem  requiring  knowledge  in  several  specialized  areas 
such  as  optics,  optics  control,  estimation  theory,  control 
theory,  etc.  This  would  be  a  very  demanding  task,  directly 
exploiting  the  tracking  filter  designed  in  this  research 
and  extended  as  discussed  above. 


Computer  Support 

The  pace  and  extent  of  this  research  effort  was  signif¬ 
icantly  retarded  by  the  computer  support  available  at  AFIT. 
An  estimated  30  percent  more  work  could  have  been  accomp¬ 
lished  with  a  more  reliable  and  responsive  system.  The 
authors  encourage  any  steps  to  make  AFIT's  computer  support 
better.  Possible  recommendations  include  more  hardware 
support  ti.e.  it  is  unbelievable  that  AFIT  has  only  one 
plotter),  more  personnel  support,  and  AFIT's  own  dedicated 
computer . 
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